#----------------------------------------------------------- #Script for stability of equilibria #with ice-albedo feedback. It also provides #a utility for finding all the equilibria in #a given temperature interval, and allows for #exploration of hysteresis phenomena. Two different #ways are given to make a hysteresis loop: integrating #the time-dependent equations, or using a (somewhat less #reliable) Newton-method solver. #----------------------------------------------------------- #Data on section of text which this script is associated with Chapter = '3.**' Figure = 'None' # #ToDo: **It would be nice to be able to take input for # the initial condition interactively, using a mouse click # # **Make L and prad into L(t) and prad(t), so that # the script can be used to illustrate hysteresis # import math from ClimateUtilities import * import phys #Define albedo and OLR functions here #Radiating pressure pRad, in mb, is passed as a #parameter. Note that the param argument could be #used to pass an object containing many values. def OLR(T,pRad): return phys.sigma * (T**4.)*(pRad/1000.)**(4.*2./7.) #Parameters of albedo function are globals def albedo(T): if T < T1: return alpha_ice elif (T >= T1)&(T<=T2): r = (T-T2)**2/(T2-T1)**2 return alpha_0 + (alpha_ice - alpha_0)*r else: return alpha_0 #Thermal mass function, dE/dT def M(T): return 50.*2000.*1000. #Rate function for use with integrator # t is not always used, but it has to be # in the argument list anyway. The # parameters needed (solar constant and # radiating pressure) are passed through the # object params. Note that this object # can also be used to pass L and pRad as # functions of time, for hysteresis studies def NetFlux(t,T,params): L = params.L pRad = params.pRad return ((L/4.)*(1.-albedo(T)) - OLR(T,pRad) )/M(T) #Function to integrate the time dependent #temperature equation and return lists of #time and temperature suitable for plotting. #The argument is the initial temperature. # #If the optional Lfunction argument is present, it #makes the solar constant L a function of time, #so that a time-dependent hysteresis loop can #be generated. tmax (days) is a global def Tseries(Tinit,Lfunction=None): # Set up the integrator dt = 24.*3600. # Time step in seconds model = integrator(NetFlux, 0.,Tinit,dt) model.setParams(params) #Do the integration tList = [] TList = [] t = 0. count = 0 while t < tmax*24*3600.: if not (Lfunction == None): params.L = Lfunction(t/(24.*3600.)) t,T = model.next() if count%50 == 0: tList.append(t/(24.*3600.)) TList.append(T) count += 1 return tList,TList #Dummy object for passing parameters. The #class Dummy is defined in ClimateUtilities params = Dummy() #Globals L = 1500. params.L = L alpha_ice = .6 alpha_0 = .28 T1 = 260. T2 = 295. pRad = 640. params.pRad = pRad tmax = 4000. # c = Curve() first = True for Tstart in [200. ,220.,240.,260.,265.,270.,280.,300.,310.,320.]: t,T = Tseries(Tstart) if first: c.addCurve(t,'Time') first = False c.addCurve(T,'Tstart=%f'%Tstart) c.Xlabel = 'time' c.Ylabel = 'Temperature' c.PlotTitle = 'Time series for various initial temperatures' w = plot(c) #------Hysteresis loop in the time-dependent system----- #This function ramps up L from an initial value up #to a maximum and then back down to a minimum, and repeates #with a given period. #Note t is in days. Period (days) is a global def Lfun(t): Lmax = 3000. Lmin = 500. return Lmin + .5*(Lmax-Lmin)*(math.sin(2.*math.pi*t/period)+1.) period = 8000. #Try longer and shorter periods. With longer #periods you get closer to the steady state for #the instantaneous value of L. You have to #go up to 64000 or more before things start to #look really close to the hysteresis diagram tmax = 3.*period #Integrate for 3 periods c1 = Curve() t,T = Tseries(200.,Lfun) #This time, we want to plot the #time series in L-T space, to #see the effect of the hysteresis #loop LList = [Lfun(tt) for tt in t] c1.addCurve(LList,'L') c1.addCurve(T,'Temperature') c1.Xlabel = 'Solar Constant (W/m**2)' c1.Ylabel = 'Temperature' c1.PlotTitle = 'Time dependence in bifurcation space' plot(c1) # #-------One way to compute a hysteresis loop ------ #Section to compute a temperature hysteresis loop in which a parameter #such as pRad is decreased and then increased again to its original #value. This implementation uses a Newton solver to get a steady #state for each value of the control parameter #ToDo: Replace or augment the following with a loop that #compiles a list of all equilibria in a given interval, and #marks them as stable or unstable. Doing the graphics to #show the stable and unstable branches for this data is a little #tricky, though. #The current implementation is written in a rather awkward #and ad-hoc way, and should be cleaned up. This calculation works #by using Newton's method to find a solution to the nonlinear #equations determining balance. The problem is that there isn't any #general way to determine which equilibrium is the right one to #switch to, so the following code has some tricks in it based #on what we know about the solutions of the system. It won't really #be reliable for a general system. # #A more reliable approach to hysteresis is to make the hysteresis #parameter time-dependant in the previous time-evolution code given #above. If the parameter is made slowly varying enough, then #the correct stable solutions will be automatically picked out. #One can also then experiment with making the time-dependence faster, #to see what happens. Note that in a hysteresis loop, the #solution can never "pass through" an unstable equilibrium. It always #goes directly from one equilibrium to a stable one. # #Perhaps this hysteresis code should be made into a separate script. # # We treat the parameters like pRad and L as function arguments rather # than globals, since that makes it easier to change their values # inside some other function. We don't need to change the # constants of the albedo function very much, so it works well # enough to leave them as globals. # #ToDo: Clean up the business of switching which control parameter #we are varying. Perhaps just allow for a general p-L variation #Reset the globals #Globals L = 1500. params.L = L alpha_ice = .6 alpha_0 = .28 T1 = 260. T2 = 295. pRad = 640. params.pRad = pRad tmax = 4000. # pStart = 1000. gList1 = [pStart-10.*i for i in range(70)] gList2 = gList1[:] #Makes a new copy of the list, which we will reverse #and append to get our list for driving the hysteresis #loop gList2.reverse() #Reverses list in place! gList = gList1 + gList2 #Our list of pRad, which decreases then increases #Do the same to make a loop of L, if we want to do a hysteresis #loop over solar constant instead Lstart = 900. LList1 = [Lstart + 10.*i for i in range(100)] LList2 = LList1[:] LList2.reverse() LList = LList1+LList2 def f(T,params): L = params.L pRad = params.pRad return (1.-albedo(T))*L/4. - OLR(T,pRad) root = newtSolve(f) root.setParams(params) pRad = gList[0] T = root(200.) Thyst = [] cHyst = Curve() #cHyst.addCurve(gList) #params.L = 1370.*.7 #for pRad in gList: # params.pRad = pRad cHyst.addCurve(LList) params.pRad = 670. for params.L in LList: T = root(T) for trial in range(1,400): if T != 'No Convergence': Thyst.append(T) if root.deriv(T,params) > 0.: print 'Warning: branch went unstable at T=%f'%T break # Get out of the loop. We're done else: #If convergence failed, that means the #solution we were tracking probably failed #to exist. Try again with different initial #guesses until we find another solution #Note that this code does not switch branches #if the branch goes unstable, and does not guarantee #that the new branch is stable. If an unstable branch #is encountered, though, a warning is issued. guess = Thyst[-1] + .5*trial*(-1)**trial T = root(guess) if T == 'No Convergence': print 'Solution cannot be continued' break #Plot hysteresis graph cHyst.addCurve(Thyst) plot(cHyst) #-------Yet another way to make the bifurcation diagram------ # #ToDo: Find a way to mark unstable branches # #This is yet another way to get hysteresis information. #This time, we find all solutions, stable and unstable, #for each value of the control parameter. This is #done by generating initial guesses by a coarse scan #of temperature space. With this method,we can #see the unstable branch explicitly. Note we could also #keep track of which branches are stable, by looking #at root.deriv(x). # #The approach used is brute-force and very inefficient, but #it's good enough for such a simple problem. For the simple #control parameter used here, the method in IceAlbedoZeroD.py #is a lot simpler, but the technique below is useful when #the dependence of the radiation balance on the control #parameter is too complicated to use that method. As an example, #you can try modifying the following code so that you do the #bifurcation diagram with ice albedo as the control parameter. #With this method, we don't need to vary the parameter up then #down. We get all solutions on the first pass. LList = [Lstart + 5.*i for i in range(200)] gList = [pStart-5.*i for i in range(140)] TList = [] #This time, it will be a list of lists, #because of multiple solutions params.pRad = 670. for params.L in LList: guesses = root.scan([200.,330.],500) solutions = [] for guess in guesses: solutions.append(root(guess)) if solutions[-1] == 'No Convergence': solutions.pop() TList.append(solutions) #Now process the solution list into a single two-column #table suitable for scatter-plotting Lmerge = [] Tmerge = [] for i in range(len(LList)): for T in TList[i]: Lmerge.append(LList[i]) Tmerge.append(T) c = Curve() c.addCurve(Lmerge,'L') c.addCurve(Tmerge,'T') c.scatter['T'] = True #Do scatter plot, don't connect with lines plot(c)