The above answers calculate equilibrium temperatures rather than the mean surface temperature. The proper spherical cow to start with applies the Stefan-Boltzmann law at each point on the surface to get the mean surface temperature. For a tidally-locked blackbody sphere (albedo = 0; emissivity = 1), this gives the following:
R code:
## Load Package to uniformly distribute lat/lon points on a sphere ##
require(geosphere)
## Steffan-Boltzmann Law ##
SBlaw <- function(I, alpha = 0.0, epsilon = 1, sigma = 5.670373e-8){
(I*(1-alpha)/(epsilon*sigma))^.25
}
## Calculate intensity of sunlight at each lat/lon ##
# The light is brightest at lat = 0, lon = 0 (max = 1362 W/m^2)
# We need to convert lat/lon to radians for R's cos function
# Irradiance cannot be negative, so a lower bound is set at zero
Imax = 1362
Npts = 1000000
LonLat = randomCoordinates(Npts)*pi/180
Irrad = pmax(0, Imax*cos(LonLat[, "lon"])*cos(LonLat[, "lat"]))
## Mean Surface Temperature ##
mean(SBlaw(Irrad))
## Equilibrium Temperature ##
SBlaw(mean(Irrad))
Results:
> ## Mean Temperature ##
> mean(SBlaw(Irrad))
[1] 157.4246
>
> ## Equilibrium Temperature ##
> SBlaw(mean(Irrad))
[1] 278.333
If you set alpha to the usual value of 0.3 (albedo), you will get ~144 k and ~255 K respectively. As the energy is smoothed over the surface of such a sphere, the mean surface temperature will approach the equilibrium temperature. The main insight hidden by using the "usual" approach is that you can get a very large changes in average temperature without putting any additional energy to the system (ie, by changing the surface distribution of energy/emissivity/albedo).
This cow is still a bit too spherical for my taste. It would be great if someone could expand on this to include rotation and surface distributions of albedo, and heat capacity. I will see about adding it later if I have time.
Edit:
Ok, I gave it a shot but really don't know how to reasonably model the energy storage in the surface for this. In case it helps anyone, here is what I could generously call a framework for a 3D rotating spherical object with latitude dependent albedo, but without atmosphere.
Code to plot the progress (can be ignored if plots is set to FALSE in the main script):
## A function to plot the progress; does not affect results ##
plotFunc <- function(Ncolors = 100, colPallet = rev(rainbow(Ncolors + 1, end = 4/6))){
if(j %% 100 == 0){
col1 = colPallet[as.numeric(cut(prev$Irrad, breaks = seq(0, 1370, length = Ncolors)))]
col2 = colPallet[as.numeric(cut(prev$Temp, breaks = seq(0, 400, length = Ncolors)))]
col3 = colPallet[as.numeric(cut(albedo, breaks = seq(0, 0.9, length = Ncolors)))]
par(mfcol = c(3,2))
plot(prev$lon, prev$lat, pch = 16, cex = .5, col = col1, panel.last = grid(),
xlab = "", ylab = "", main = "Insolation (W/m^2)")
map(plot = T, fill = F, add = T)
image.plot(matrix(rnorm(10)), breaks = seq(0, 1370, length = Ncolors+2),
col = colPallet, legend.only=T, horizontal=T)
plot(prev$lon, prev$lat, pch = 16, cex = .5, col = col3, panel.last = grid(),
xlab = "", ylab = "", main = "Albedo (% Reflected)")
map(plot = T, fill = F, add = T)
image.plot(matrix(rnorm(10)), breaks = seq(0, 90, length = Ncolors+2),
col = colPallet, legend.only=T, horizontal=T)
plot(prev$lon, prev$lat, pch = 16, cex = .5, col = col2, panel.last = grid(),
xlab = "", ylab = "", main = "Temperature (K)")
map(plot = T, fill = F, add = T)
image.plot(matrix(rnorm(10)), breaks = seq(0, 400, length = Ncolors+2),
col = colPallet, legend.only=T, horizontal=T)
plot(colMeans(tempHistory[,1:cnt]), type = "l", xlab = "Time Step",
main = "Mean Surface Temperature", ylab = "Temperature (K)", lwd=2)
dens = density(TempSurr)
hist(prev$Temp, freq = F, col = "Grey", xlab = "Surface Temperature (K)",
main = "Distribution of Surface Temperatures",
breaks = seq(0, max(TempSurr, prev$Temp), length = 40))
lines(dens, col = "Red", lwd=3)
abline(v = c(mean(TempSurr), mean(prev$Temp)), col = c("Red", "Black"),
lwd =3, lty = c (1,2))
}
msg = cbind(dT = c(range(dT), mean(dT)),
Temp = c(range(prev$Temp), mean(prev$Temp)),
TempSurr = c(range(TempSurr), mean(TempSurr)))
rownames(msg) = c("min", "max", "mean")
print(paste("Day = ", d, " Solar Angle = ", j))
print(msg)
}
The simulation code:
## Load Packages ##
require(geosphere)
require(maps)
require(fields)
## Choose whether to make the plots ##
plots = TRUE
## Steffan-Boltzmann Law ##
SBlaw <- function(I, alpha = 0.0, epsilon = 1, sigma = 5.670373e-8){
(I*(1-alpha)/(epsilon*sigma))^.25
}
## Initialize misc parameters ##
# The coordinates should be spread uniformly over the sphere
# The light will be brightest at lat = 0, lon = 0 (max = 1362 W/m^2)
# The object will rotate relative to sun at w ~ 0.004 degrees lon per sec
# Use a simple albedo model that is a function of latitude
# c is a "thermal resistance" constant. Temp can only rise by c*(Radiation Temp - Current Temp)
Imax = 1362
w = 7.2921150e-5*180/pi
LonLat = as.data.frame(regularCoordinates(50))
c = .01
# S-B law parameters
epsilon = 1
sigma = 5.670373e-8
albedo = abs(LonLat$lat/100)
#albedo = albedo[order(abs(LonLat$lat))]
# The model will update once every x*w seconds for nDays
tStep = 5*60
nDays = 5
offsets = seq(0, 360, by = tStep*w)
prev = cbind(LonLat, Irrad = 0, Temp = 0)
tempHistory = matrix(nrow = nrow(prev), ncol = nDays*length(offsets))
cnt = 0
for(d in 1:nDays){
for(j in 1:length(offsets)){
# We need to convert lat/lon to radians for R's cos function
# Irradiance cannot be negative, so a lower bound is set at zero
IrradIn = pmax(0, Imax*cos((LonLat$lon + offsets[j])*pi/180)*cos(LonLat$lat*pi/180))
IrradOut = epsilon*sigma*prev$Temp^4
IrradNet = (1- albedo)*IrradIn - IrradOut
TempSurr = SBlaw(pmax(0, IrradNet))
# The actual change in temp is a function of the imbalance between the current
# temp and that it should be at if at equilibrium with the incoming radiation.
# This most likely means nothing, it is a placeholder!!!
dT = c*(TempSurr - prev$Temp)
# Update Temperatures + Irradiation
prev$Temp = prev$Temp + dT
prev$Irrad = IrradIn
# Store the temperatures
cnt = cnt + 1
tempHistory[, cnt] = prev$Temp
if(plots){ plotFunc() }
}
}
If anyone has any ideas for modeling the storage of energy in the surface of this object in a simple way, please share.
You can see my attempt gave interesting results. The average temperature actually did not change from the "tidally-locked" object, but the distribution did. This can be seen in the bottom plot (red ~tidally locked distribution; histogram = current model).
