|Geo 161 Lab 4: Analyzing model output
Objectives: Paleoclimatology often uses climate models for interpretation of geochemical proxy data. In class, we’ve talked about using climate models to understand how Earth goes into and comes out of a ‘snowball’ state, to understand the maintenance of high latitude warmth and suppression of latitudinal temperature gradients during ‘greenhouse’ climates, and to explain how small changes in Northern hemisphere insolation might ultimately trigger the growth and decay of continental ice sheets (to name just a few). Sometimes such models can seem a bit like a black box -- both in terms of how the model actually works and in terms how one might go about viewing or analyzing model output. For this lab, my goal is to familiarize you with climate model output using the intermediate complexity model cGENIE (Grid ENabled Earth System Model) through both MATLAB and a software program called Panoply.
Background: Some of this model output will be time series data similar to the types of datasets you’ve already analyzed, which we can easily load and plot in MATLAB. However, we will also look at spatially gridded data (remember, 3D models break up the area of the Earth into individual grid cells and represent average conditions within each grid). Model output files of gridded data are saved for this model in the netCDF file format. NetCDF stands for Network Common Data Format, which is a self-describing machine-independent data format for storing array-oriented data (for ‘array-oriented’, just think of data in rows and columns, like an mxn matrix in MATLAB). A netCDF file not only contains the raw numeric data, but also identifying information about that data. We’ll take a look at netCDF files in MATLAB briefly, but to make quick plots and nicely visualize our results without a lot of code, we’ll use a software program called Panoply that has been installed on all the Desktop computers in Watkins 2101.
Note: If you are using your own computer, you can download Panoply (it’s free!) for your operating system. (Go to http://www.giss.nasa.gov/tools/panoply/ and simply make sure that you also have a Java SE 7 or 8 runtime environment installed.
We will take a look at two different sets of experiments that I’ve run using the cGENIE model. I’ve uploaded model results on ilearn.ucr.edu under Assignments where each folder represents an experiment and contains a subfolder named biogem, which then contains individual results files. The files with .res filename endings are text files that we can load easily into MATLAB using the same commands we’ve used previously. The .nc filename endings are the netCDF file type for spatial model output. First, create a folder in your MATLAB directory called cgenie_output. (Remember, you can do this using the following command)
Next, download all the experiment folders and add them to the cgenie_output directory within your MATLAB working directory.
4.1 Ocean Circulation in cGENIE
Changes in deep ocean circulation through time can influence global heat transport, air-sea exchange of CO2 and other gases, and biogeochemical cycles. As we move into discussion of the glacial-interglacial cycles of the Pliocene and Pleistocene (and especially the asymmetrical, ‘saw-tooth’ shaped pattern of slow cooling and rapid warming that characterized the late Pleistocene) you will notice an increased focus on the behavior of the global overturning circulation (also called the meridional overturning circulation or thermohaline circulation). Correspondingly, lots of attempts have been made to model the overturning circulation and test its stability in response to changes in global climate.
The cGENIE model includes a three-dimensional ocean model, so we can evaluate patterns and strength of the overturning circulation. There is also an ocean biogeochemical module (called ‘BIOGEM’) coupled to the physical circulation model, through which we can track the spatial and temporal concentration of a number of different tracers – various dissolved species, like dissolved inorganic carbon (DIC) or characteristics of water masses like temperature and salinity. Through BIOGEM we can also apply ‘forcings’ to the system. These are time and/or spatially dependent fluxes or restored-to boundary conditions that can be prescribed for the various tracers. In this first exercise we will look at two kinds of forcings – one of a colored dye tracer and another of freshwater (actually a negative salinity flux).
In cGENIE, locations are specified on a grid of (i,j,k) coordinates. The i and j values are equivalent to longitude and latitude, respectively, and the k values are equivalent to depth. These experiments are run on a 36x36x16 grid, meaning that the longitude and latitude resolution is 1° while the ocean is divided into 16 depth levels. Depth levels cover narrower ranges at the surface and increase in size with depth – the surface ocean is an 80 m thick layer, while the deepest ocean layer is >1000 m thick. The continental configuration and ocean bathymetry can be altered for paleo-experiments, but today we will look at experiments using a modern continental configuration only.
In order to trace the deep ocean circulation in cGENIE, the first experiment we will analyze is called LAB_2.colorinjection and has had a pulse of red dye applied at a single point in the North Atlantic ocean south of Greenland at intermediate depths over the course of the experiment. This pulse of dye is called a flux forcing in units of mass per year.
In this experiment, the model has been run for 100 years from a restart. The restart means that the model was previously run for 10,000 years to stabilize ocean circulation and atmospheric CO2. This is a ‘closed’ system experiment for all considered tracers, which just means that tracers may be advected (moved about) within the ocean, but their total concentrations remain unchanged (no addition or removal).
4.1.1 Viewing netCDF results.
You can see how the red tracer distribution evolves (and thus visualize the ocean circulation) by looking at the spatial data contained in the netCDF files fields_biogem_2d.nc and fields_biogem_3d.nc, for the two- and three-dimensional datasets, respectively. Actually, these are spatial datasets for a number of different variables that evolve with time. In MATLAB,
Will display the contents of the netCDF file fields_biogem_3d.nc. Scroll to the top of the list that appears in the Command Window and you’ll see some basic attributes about the file, including the dimensions for all the variables. Notice the 36x36x16 spatial dimensions, and also notice the time dimension = 7. In other words, there are 7 time-slices saved over the course of this 100-year experiment.
Scroll down through the list of variables to see what else has been saved for this experiment. Note the variable ocn_colr. This is our red dye tracer, and you can see that it is saved as a four dimensional matrix (36x36x16x7) or longitude x latitude x depth x time.
We’ll come back to look at netCDF data in MATLAB later, but first let’s use a quicker way to visualize spatial results. In your Current Folder window, double click on the folder cgenie_output, then LAB_2.colorinjection, then biogem, and finally right-click on the file fields_biogem_3d.nc. Select the option to ‘Open Outside MATLAB.’ This should launch the program Panoply. (If it doesn’t, check that Panoply is installed on the machine. You can also try to open Panoply first and then choose File>Open File from the dropdown menu and navigate to where LAB_2.colorinjection/biogem/fields_biogem_3d.nc is saved in your MATLAB/cgenie_output folder).
Look at results for the variable ocn_colr by simply clicking on the variable name on the left-hand side of the Panoply window and selecting Create Plot. You will get an option to create different types of plots. For now, try the default Longitude-Latitude plot and select Create.
By default, Panoply will overlay the grid with a modern map image. In cGENIE, grid cells for land are colored grey. Notice the difference in detail between the map outline and the cGENIE continental grid.
Question 4.1: Do you think that this model would do a good job of representing the Indonesian Throughflow? Why or why not? Recall, we discussed this ocean gateway when considering causes of the Pliocene transition and onset of Northern Hemisphere glaciation. Hint: look north of Australia.
Un-check the box labeled Interpolate in order to get a better idea of the resolution in cGENIE. When you create a new map in Panoply, it will automatically start with Year 1 and the surface ocean (Z-level mid depth), which it will describe both as an integer layer number and a depth mid-point (in meters). NOTE: Panoply counts the layers differently than the model; it is calling the surface layer 1, while in the model, the surface layer is actually indicated by k=16. Click through the ocean layers and experiment years to watch the movement of the dye tracer.
Question 4.2: Do you see any evidence of the red dye tracer in the central Pacific Ocean during this experiment? Does this surprise you? Why or why not?
You can also view the red dye tracer as a water column profile. Back in the main Panoply window, double-click the ocn_colr variable again, but this time, create a Latitude-Vertical plot. By default, what you’re looking at is a global average across all longitudes. To focus on the Atlantic, un-check the box Avg and for Longitude of the t-grid, go to a value of 22 out of 36 (this is the value specifying the longitude where the dye was first injected). Now you can similarly move through the years of the experiment to watch the dye propagate.
Question 4.3: Does the dye ‘spread out’ more in the surface ocean or the deep ocean? What do you think causes this pattern?
In addition to using the dye injection to look at ocean circulation in the model, you can also plot the overturning circulation directly from the 2D netCDF output using the variable phys_opsi, which is the global overturning streamfunction, or phys_opsi_a, which is the overturning streamfunction in just the Atlantic. An overturning streamfunction tells you the volume transport of water using the unit Sverdrups (Sv), which is equivalent to transport of a million cubic meters per second. These output variables should help provide a visualization of the large-scale ocean circulation that is driving the movement of the dye tracer. The positive and negative just tell you the direction (clockwise or counterclockwise, respectively) of the overturning circulation, while the absolute value of the number tells you the magnitude of the volume transport.
I have also provided you with an experiment that modifies the location of the red dye input called LAB_2.colorinjection.SAtl, which adds the same magnitude of dye injection into a point at the same depth and same longitude but in the South Atlantic (i,j,k)=(22,4,8). Try visualizing the output of this experiment compared to the previous experiment (hint: open another Panoply window of the same tracer and compare them side by side).
Question 4.4: Does the dye move similarly over similar timescales from its point of origin (spatially or with depth)? What physical processes could explain the different paths traced by the dye?
4.1.2 Sensitivity of overturning circulation to freshwater forcing
Instead of adding a dye tracer, which is just being passively advected by the ocean circulation, it is possible to add fresh water to assess the sensitivity of the Atlantic Meridional Overturning Circulation (AMOC) to collapse. As we will discuss in the next few weeks, major changes in AMOC have been invoked to explain glacial-interglacial climate changes as well as shorter-term climate events, particularly during deglaciation. This is a classic test for ocean models, called a ‘hosing’ experiment.
In the experiment LAB_2.hosing, the location of the freshwater forcing is in the North Atlantic at the same site as injection of the red dye in the previous experiment. However, this time the flux is at the surface instead of intermediate depths. The idea is that we can represent an instance of glacial meltwater injection to the surface of the North Atlantic Ocean. Because cGENIE has a fixed volume constraint (the ocean is not gaining or losing volume), this freshwater flux is actually achieved as a negative salinity flux (think of it as removing salt from seawater). I won’t go through the conversion, but the negative salinity flux applied is equal to a sustained 0.09 Sv freshwater input. In general, fluxes of 0.05 Sv and larger (and to specific places) are applied to models in order to cause the AMOC to collapse.
In order to investigate the size of the freshwater flux that can be sustained without collapsing the AMOC, I have provided you with an additional set of experiments.
These, respectively, double, half, and quarter the magnitude of the freshwater flux forcing.
View results of these experiments by opening the file fields_biogem_2d.nc from each experiment in Panoply and making plots of the Atlantic Overturning streamfunction (phys_opsia). Click through the time-slices and pay attention to the North Atlantic overturning strength.
Question 4.5: What is happening to the North Atlantic overturning strength as a consequence of freshwater forcing in the first experiment LAB_2.hosing? Do any of these experiments cause an apparent collapse in the AMOC? Which one(s)?
We can also look at the overturning strength as time-series data, which we can easily plot in MATLAB. Go back to the MATLAB Command Window and type
This dataset has 5 columns: Year, global min overturning (Sv), global max overturning (Sv), Atlantic min overturning (Sv), and Atlantic max overturning (Sv). The min and max just refer to the direction (counterclockwise or clockwise overturning). To plot the North Atlantic overturning streamfunction, use the 5th column.
Load the overturning streamfunction time series from the other hosing experiments, and remember to give them each a different name (i.e. opsi_hosingx0p25, etc.)
Question 4.6: Make a plot that compares the maximum Atlantic overturning in each experiment.
There are important climatic implications of disrupting the AMOC, which is why changes in AMOC are often invoked to explain decadal to centennial climate change, especially at the high latitudes. Let’s try visualizing the importance of the AMOC for surface air temperatures, which should give you some indication of the importance of the AMOC in transporting heat to the North Atlantic region. To more easily assess what’s going on during these experiments it is possible to create an anomaly (difference) map in Panoply:
First, create a plot of atm_temp (surface air temperature) from the 2D netCDF file using one of the hosing experiments (pick one that you think is likely to have noticeable effects on climate). Next, with atm_temp still selected (and the first plot window still open), click on the ‘Combine Plot’ icon. A dialogue box will appear and ask you to select a plot to combine with the new one. Make sure that the name of your first plot window is selected/highlighted and then click ‘Combine.’
You now have a plot window that by default shows the difference between two identical time slices. The two different slices are labeled Array 1 and Array 2. Select to plot a difference map of Array 2 – Array 1 from the drop down menu. Then, keep Array 1 fixed to the initial year (Year 1 = centered on 0.5 years) and vary the year in the second array (Array 2). To really see the difference, you’ll need to go to the Scale tab and select Fit to Data or Center on Zero. You can also play around with the colors used under the drop down menu for Color Table and also invert colors to get a plot that highlights the pattern of changing surface temperature.
Question 4.7: What is the consequence of ‘hosing’ the North Atlantic for surface air temperatures? Do surface air temperatures change the same everywhere?
4.2 Fossil Fuel CO2 release and ocean acidification
Much of the motivation for studying paleoclimate comes from the desire to provide context for our modern fossil fuel CO2 emissions ‘experiment.’ In these next experiments, we will consider some of the impacts on the marine carbon cycle of anthropogenic CO2 release.
Note that cGENIE is not the model to use if you’re interested in constraining climate sensitivity (or how much warming will result per doubling of atmospheric CO2). Rather, cGENIE has a default climate sensitivity of 2.64°C warming per CO2 doubling. However, cGENIE has a detailed representation of marine carbonate chemistry, so this is an ideal model for thinking about changes in marine CO2 uptake and the consequences in terms of ocean chemistry.
In the experiments you’ll be looking at, the rate of calcification by plankton living in the surface ocean is responsive to ocean acidification and saturation state (i.e. it takes into account CO2-calcification feedbacks which will additionally interact with climate). Remember, in class we discussed that decreased carbonate production increases ocean CO2 uptake, which acts as a (possibly counterintuitive) negative feedback for atmospheric CO2.
4.2.1 Time-series plotting in MATLAB
The experiment LAB_3.CO2emissions.1yr involves a release of CO2 to the atmosphere (globally uniform) of 1000 Pg C over a single year. Remember, 1 Pg = 1 billion tons of carbon. The isotopic composition of that carbon (its δ13C) is -27‰ to represent fossil fuel CO2.
First, we will look at time series results in response to this CO2 forcing in MATLAB. Double click on the experiment name and then the biogem subfolder to view all the .res files included for this experiment. First, we will load and plot the evolution of atmospheric CO2.
This will plot the 3rd column in the pCO2.res file, which has units of μatm (hence the multiplication by 1e6 above to get values that look more familiar). Notice that the extremely rapid input of CO2 over just one year, combined with the saving of average values for each year, means that the first year atmospheric CO2 concentration is already well above today’s value although this experiment was run from pre-industrial CO2 of 280 ppm!
Try looking at changes in surface ocean pH compared to atmospheric pCO2. First, load the pH data:
Next, we can plot two different variables on two different y-axis using the plotyy function – useful for making comparisons between the shape of each plot.
The above will generate a plot with two y-axis, with atmospheric CO2 on the first, left-hand axis, and surface ocean pH on the second, right-hand axis. The output to the plotyy function, [Ax,H1,H2] returns something called axes handles, which we can use to modify the individual axes. For instance, try reversing the direction of the y-axis for surface ocean pH using the set command.
Now you can compare the evolution of atmospheric CO2 and surface ocean pH in response to the 1 yr CO2 pulse.
Question 4.8: Which variable shows a faster recovery? Why do you think this might be the case? (Hint: where is the atmospheric CO2 going)?
I have also provided you an experiment with a slower rate of carbon injection LAB_3.CO2emissions.100yr. Now, instead of 1000 Pg C applied in a single year, emissions are at a rate of 10 Pg C per year over 100 years. In other words, the total emissions amount is the same in both experiments. Load the CO2 dataset and pH dataset from this experiment and make another plot equivalent to the one you made for the 1 yr emissions case. (Reminder: if you use the figure command before making the next plot, you will not override your previous plot).
Question 4.9: What are the major differences between the 1 yr and 100 yr release in terms of these two variables?
In addition to simple pulse forcings of CO2, we can force the model with historical and future emissions scenarios. The experiment LAB_3.historical prescribes the CO2 emissions to match the observed historical concentration profile from year 1765 until 2010 while the experiment LAB_3.future continues where the historical experiment leaves off but uses an Intergovernmental Panel on Climate Change emissions scenario out to the Year 2100, with total emissions of 2180 Pg C.
Load the CO2 and pH datasets for each of these experiments and then combine them into single datasets.
Plot these data using plotyy.
Question 4.10: Do you think the more realistic historical-future scenario more closely resembles the 1 yr release or 100 year release experiment? What might be some consequences of using a 1 yr pulse release experiment to simulate the impacts of modern CO2 release on ocean carbon chemistry?
4.2.2 NetCDF results in MATLAB
You can look at spatial impacts of CO2 release on things like surface ocean pH using the netCDF results from these experiments in Panoply (and I encourage you to do so if you have time). You can also generate similar spatial plots in MATLAB, but we won’t do that here since it involves a lot more code. Instead, another useful aspect of analyzing spatial data in MATLAB is that you can analyze or create time-series records for individual grid cells.
From the ilearn.ucr.edu page, download the m-file called load_netcdf_basic.m to your MATLAB working directory (NOT within the cgenie_output folder). Double-click on the file name in the Current Folder window to view the file in the Editor window. This is just a short function that will allow you to load variables from either the 2D or 3D cGENIE netCDF files. The format for the function is [output]=load_netcdf_basic(EXP,VAR,DIM) where output is the output matrix (name it whatever you like) and the required inputs are EXP for the experiment name, VAR for the variable name and DIM for the dimension. Note that the inputs for EXP and VAR are strings (in quotes) and DIM is an integer (either 2 or 3).
Remember, the ncdisp command will tell you what variables are contained in each file. You can also load these files into Panoply to see what variables are available for plotting.
Will give you a four-dimensional array of pH data from the experiment LAB_3.future. The first two dimensions are longitude and latitude, the third is ocean depth, and the fourth is time. In other words, you have an array containing pH values for the full 3D grid for each of the time slices saved (there are four in this experiment at years 2019.5, 2049.5, 2099.5, and 2159.5).
Another important variable to assess is the surface ocean saturation state with respect to the mineral aragonite. Remember, aragonite is calcium carbonate, but it’s a form of the mineral used by coral to make their skeletons. Aragonite is physically stronger than calcite, but it is more susceptible to dissolution. Ocean surface waters in which aragonite becomes under-saturated (Ω < 1.0) are regarded as a critical threshold for organisms making aragonite shells and skeletons. Load the surface ocean aragonite saturation dataset and also the years for which we have time-slice data:
Let’s pick one site near the equator, try (15,15), and plot the trajectory of the aragonite saturation state at that site.
The permute command simply switches the order of the dimensions so that we can plot against time.
Try looking at a few more individual sites (you choose which ones!)
Finally, we can also compare these individual sites to the average value saved for the entire surface ocean in the time series .res files.
You can also adjust values relative to zero if you want to look at the size of the change rather than the absolute values.
>> hold on
Question 4.11: If you pick a high and low latitude site, which sees a relatively larger drop in aragonite saturation state?