diff --git a/Time-Stepping Naked Planet Model/Grading/GWmod1.py b/Time-Stepping Naked Planet Model/Grading/GWmod1.py deleted file mode 100644 index 90356e0..0000000 --- a/Time-Stepping Naked Planet Model/Grading/GWmod1.py +++ /dev/null @@ -1,76 +0,0 @@ -#My first tijme making a plot for climate modeling. -#Hopefully, this works. -#Note that using spreadsheets for this course is optional; -#the only way to get credit is with Python3! - -"""The background for this material can be found in Sections 2 and 3 of Part I -of this class. Joules are units of energy; Watts are energy flow (J/s). -The temperature of a planet is determined by balancing energy fluxes in and out of a planet. -Incoming solar heat is determined by L*(1-albedo)/4, -and outgoing infrared is calculated as epsilon * sigma * T^4.""" - -"""The heat capacity of the surface planet depends on the water thickness. -Using the fact that 1 gram of water warms by 1 K with a 4.2J (1 cal) of heat, -and the density of water, to calculate how many Joules of energy it takes -to raise the temperature of 1 m2 of the surface, given some water depth -in meters.""" - -"""The goal is to numerically simulate how the planetary temperature of a -naked planet would change through time as it approaches equilibrium -(the state at which it stops changing, which we calculated before). -The planet starts with some initial temperature. -The “heat capacity” (units of Joules / m2 K) of the planet is set by -a layer of water which absorbs heat and changes its temperature. -If the layer is very thick, it takes a lot more heat (Joules) to change the temperature. -The differential equation you are going to solve is:""" -#d(HeatContent)/dt = L*(1-alpha)/4 - epsilon * sigma * T^4 -"where the heat content is related to the temperature by the heat capacity:" -#T[K] = HeatContent [J/m2] / HeatCapacity [J/m2 K] -"""The numerical method is to take time steps, extrapolating the heat content -from one step to the next using the incoming and outgoing heat fluxes, -same as you would balance a bank account by adding all the income and subtracting all the expenditures -over some time interval like a month. -The heat content "jumps" from the value at the beginning of the time step, -to the value at the end, by following the equation:""" -#HeatContent(t+1) = HeatContent(t) + dHeatContent/dT * TimeStep -"""This scheme only works if the time step is short enough so that -nothing too huge happens through the course of the time step. -Set the model up with the following constants:""" -timeStep = 40 # years; convert to seconds. -time = 0 # years; this shall denote how much time has passed in total. -T = 0 # Kelvin, according to question demands. -waterDepth = 2000 # meters; multiply by cH2O and water density. -L = 1350 # W m-2, same as J s-1 m-2 (Solar Constant) -> multiply by no. of seconds in each timestep -albedo = 0.3 # alpha -epsilon = 1 -sigma = 5.67E-8 # W m-2 K-4 -cH2O = 4200 # specific heat capacity of water (units of 4200 J kg-1 K-1). -rhoH2O = 1000 # Density is 1000 kg m-3. -energy_per_m2 = waterDepth * cH2O * rhoH2O #in J m-2 K-1 -#heat influx and outflux must be calculated along the way. -import numpy -import matplotlib.pyplot -#Getting our lists needed with the very first variables already in: -nSteps = input("Please use 20 or more") -#apparently the autograder that the MOOC uses also mistakes strings for operators. -nSteps = int(nSteps) -time_list = [time] -temperature_list = [T] -HeatContent = 0 #only true when T = 0 Kelvin, this is the heat content of the Earth. -difference = 0 -outflux = 0 -for x in range(nSteps): - time = time + timeStep - time_list.append(time) - difference = (L * (1 - albedo) / 4 - (epsilon * sigma * T**4)) #W m-2 - HeatContent = HeatContent + difference * 365 * 24 * 60 * 60 * timeStep #J m-2 - T = (HeatContent / energy_per_m2) #Kelvin - outflux = epsilon * sigma * T**4 # as per instructions - temperature_list.append(T) -#derive answers and end with: -matplotlib.pyplot.plot(time_list, temperature_list) -matplotlib.pyplot.show() -#Give the last temperature derived by the model: -#print(temperature_list[-1], outflux) -#The tolerances in the grader for temperatures are 2-3 degrees C, -#and for the heat fluxes, +- 0.1 W/m2. diff --git a/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py b/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py deleted file mode 100644 index 8a63aea..0000000 --- a/Time-Stepping Naked Planet Model/Grading/Naked_Earth.py +++ /dev/null @@ -1,32 +0,0 @@ -import matplotlib.pyplot as plt - -timeStep = 10 # years -waterDepth = 4000 # meters -L = 1350 # Watts/m2 -albedo = 0.3 -epsilon = 1 -sigma = 5.67E-8 # W/m2 K4 -#nSteps = int(input("")) -nSteps=200 - -heat_cap=waterDepth*4.2E6 #J/Kg/deg -time_list=[0] -temperature_list=[400] -heat_content=heat_cap*temperature_list[0] -heat_in=L*(1-albedo)/4 -heat_out=0 -secs_in_yr=3.14E7 - -for t in range(0,nSteps): - time_list.append(timeStep+time_list[-1]) - heat_out=epsilon*sigma*pow(temperature_list[-1],4) - heat_content=heat_content+(heat_in-heat_out)*timeStep*secs_in_yr - temperature_list.append(heat_content/heat_cap) - heat_out = epsilon * sigma * pow(temperature_list[-1], 4) -print(temperature_list[-1]) -#print(temperature_list[-1], heat_out) -plt.plot(time_list, temperature_list, marker='.', color="k") -plt.xlabel("Time (years)") -plt.ylabel("Temperature (degrees Kelvin)") -plt.show() - diff --git a/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py b/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py deleted file mode 100644 index 4d567f2..0000000 --- a/Time-Stepping Naked Planet Model/Grading/naked_planet_tosubmit.py +++ /dev/null @@ -1,68 +0,0 @@ -# %% - - -# %% -import numpy as np -import matplotlib.pyplot as plt - -# Set up constants as given: -timeStep_years = 100 # years -waterDepth = 4000 # meters -L = 1350 # W/m² -albedo = 0.3 -epsilon = 1 -sigma = 5.67e-8 # W/m² K^4 - -# The incoming solar flux is constant: -incoming_flux = L * (1 - albedo) / 4 # W/m² - -# Calculate the heat capacity (J/m² K) -# (water density = 1000 kg/m³, specific heat = 4184 J/(kg·K)) -heat_capacity = waterDepth * 1000 * 4184 - -# Read the number of time steps from stdin -nSteps = int(input("")) - -# Set the initial conditions. -# (Initial temperature is chosen in Kelvins; 0 K is allowed though not physical.) -initial_temperature = 0.0 # Kelvin -heatContent = initial_temperature * heat_capacity # in J/m² - -# Prepare lists for time and temperature (for plotting, if desired) -time_list = [0] -temperature_list = [initial_temperature] - -# Convert the time step to seconds. -dt = timeStep_years * 365 * 24 * 3600 - -# Integration loop: update heat content and temperature each step. -for i in range(nSteps): - # Current temperature in Kelvin: - T = heatContent / heat_capacity - - # Differential equation: dHeat/dt = incoming_flux - ε·σ·(T)⁴ - dH_dt = incoming_flux - epsilon * sigma * (T**4) - - # Update the heat content over the time step: - heatContent = heatContent + dH_dt * dt - - # Append new time and temperature to lists: - time_list.append((i+1) * timeStep_years) - temperature_list.append(heatContent / heat_capacity) - -# Compute the final temperature (in Kelvin) and the corresponding outgoing heat flux. -final_temperature = heatContent / heat_capacity -final_heat_flux = epsilon * sigma * (final_temperature**4) - -# Output only the final temperature and heat flux. -print(final_temperature, final_heat_flux) - -# Plotting code -plt.plot(time_list, temperature_list) -plt.show() - - -# %% - - - diff --git a/iterative_runaway_ice_albedo_feedback.py b/iterative_runaway_ice_albedo_feedback.py new file mode 100644 index 0000000..4eac267 --- /dev/null +++ b/iterative_runaway_ice_albedo_feedback.py @@ -0,0 +1,80 @@ +import numpy as np +import matplotlib.pyplot as plt + + +### Properties ### +steps = 100 + +# Data +T_data = np.array([265, 255, 245, 235, 225, 215]) +il_data = np.array([75, 60, 45, 30, 15, 0]) +albedo_data = np.array([0.15, 0.25, 0.35, 0.45, 0.55, 0.65]) + +# Ice latitude +m1, b1 = np.polyfit(T_data, il_data, 1) + +def calc_ice_latitude(T): + return m1 * T + b1 + +# Planetary albedo +m2, b2 = np.polyfit(T_data, albedo_data, 1) + +def calc_albedo(T): + return np.clip(m2 * T + b2, 0, 1) + +### Constants and Data ### +epsilon = 1 +sigma = 5.67E-8 # W/m2 K4 + +### Energy balance equation ### +def calc_temperature(albedo, L): + return pow((L * (1 - albedo)) / (4 * epsilon * sigma), 1/4) + +### Iteration step ### +def iterate(L, initial_albedo, steps): + T = calc_temperature(initial_albedo, L) + for i in range(steps): + albedo = calc_albedo(T) + new_T = calc_temperature(albedo, L) + T = new_T + return T, albedo + +### Loop ### +l_range_bounds = [ 1200, 1600 ] + + +current_albedo = 0.15 +temp_down = [] +l_range_down = [] +L = l_range_bounds[1] +while L > l_range_bounds[0]-1: + T, current_albedo = iterate(L, current_albedo, steps) + temp_down.append(T) + L = L - 10 + l_range_down.append(L) + +current_albedo = 0.65 +temp_up = [] +l_range_up = [] +L = l_range_bounds[0] +while L < l_range_bounds[1]+1: + T, current_albedo = iterate(L, current_albedo, steps) + temp_up.append(T) + L = L + 10 + l_range_up.append(L) + + +### Input ### +# L, albedo, nIters = input("").split() +# L, albedo, nIters = [ float(L), float(albedo), int(nIters) ] + +# T, albedo = iterate(L, albedo, nIters) +# print(T, " ", albedo) + +### Plot hysteresis ### +plt.plot(l_range_down, temp_down, label="Cooling") +plt.plot(l_range_up, temp_up, label="Warming") +plt.xlabel("Solar Constant (L)") +plt.ylabel("Equilibrium Temperature (K)") +plt.legend() +plt.show() \ No newline at end of file diff --git a/Time-Stepping Naked Planet Model/naked_planet.py b/time_stepping_naked_planet.py similarity index 85% rename from Time-Stepping Naked Planet Model/naked_planet.py rename to time_stepping_naked_planet.py index 66cce23..fab0206 100644 --- a/Time-Stepping Naked Planet Model/naked_planet.py +++ b/time_stepping_naked_planet.py @@ -3,10 +3,10 @@ import matplotlib.pyplot # Inputs -timeDuration = 2000 # years -timeStep = 10 # years -waterDepth = 4000 # meters -initialTemp = 400 +timeDuration = 2000 # years +timeStep = 10 # years +waterDepth = 4000 # meters +initialTemp = 400 # K # Constants secondToYear = 60*60*24*365 @@ -41,11 +41,13 @@ for i in range(nSteps): # Save timestep for plotting time.append(time[-1] + timeStep) -print(T[-1]) +print("final temperature: ", T[-1]) +print("Initial heat out: ", heatOut[0]) +print("Final heat out: ", heatOut[-1]) matplotlib.pyplot.plot(time, T) matplotlib.pyplot.title("Time-dependent temperature of a naked planet") matplotlib.pyplot.xlabel('Time (years)') matplotlib.pyplot.ylabel('Temperature (K)') -matplotlib.pyplot.show() \ No newline at end of file +matplotlib.pyplot.show()