Thursday 23 April 2015

What is PIOMAS gice?

 In this post I'm going to explain a PIOMAS metric that I will be using from now on. I will also give a first look at this metric. The gridded data for the 'sub-grid thickness distribution' is contained in gice files, and I've finally started to use them.

I often use PIOMAS grid box effective thickness, and it is used to make plots like this, from my March 2015 status post.

PIOMAS works on a grid of curvilinear coordinates, these curves can be seen in the way lines seem to bend in the green landmass of the above image. Despite this rather daunting looking geometry, PIOMAS gridded data is provided in arrays of 360 columns and 120 rows, this provides 43200 datapoints for each month or day, depending on the source files used. Note how in the image above the 'pole hole' is not centred over geographical north. The pole hole is centred over Greenland and means the physics calculations don't have to contend with the tricky business of moving through the pole.

So those plotted points are the centre of the grid boxes, each data point marks a grid box with a specified latitude and longitude and a specified area. Grid box effective thickness is held in files called 'heff', files (h is thickness, so heff is h effective). Taking a certain grid box I can obtain its thickness by loading an array from a heff file and simply accessing what the value is, for the grid cell in question the grid box effective thickness is 1.992985m. In another array the grid box area is 580.7655km^2, using these two figures I can calculate the volume of ice in the grid box,

580.7655 * 0.001992985 = 1.15745693km^3.

If I do this for every grid box and sum the result I can work out total volume. As I know the latitude and longitude of the grid boxes I can work out which grid box falls into which region, hence the regional volume data I produce.

OK so far, but what is grid box effective thickness and what does a thickness like 1.993m actually mean?

In "A Thickness and Enthalpy Distribution Sea-Ice Model" Zhang and Rothrock explain the workings of the model that would eventually become PIOMAS. Their figure 2 shows the sub-grid thickness distribution.

Basically the model works with grid boxes, growing/melting ice, allowing ice movement compaction/thinning, but it doesn't assume that the thickness of a grid box is just 1.993m, or whatever thickness. It actually contains a thickness distribution which is contained in the gice files. For a while now Dr Zhang has made these files available, but it has only been recently with a desire to look at the multi year ice export into Beaufort and Chukchi that I have been motivated to look at that data.

Here is how the sub grid thickness distribution works. With a grid cell there is a representation of how much ice is at each thickness outlined above, in 12 bands of thickness. PIOMAS is made up of a set of equations, here's one just as an example:

No need to worry about understanding it, although this happens to be the equation that handles ridging of ice in the model as simulated movement crushes ice together creating pressure ridges.

So using the sub grid thickness distribution the above equation acts on the entire thickness distribution, in other words it is calculated at least 12 times for each thickness (the evaluation of most equations seems to be iterative) and the result is partitioned into a new sub grid thickness distribution. This means that using sub-grid parameterisation PIOMAS is able to model processes happening at smaller scales than the grid box. It is the combination of that technique and assimilation of real world data that makes the PIOMAS data such a skillful reanalysis product for sea ice.

Anyway, that grounding covered, let's take a look at some gice and heff data. Both sets of data are obtained from Dr Zhang's webpages, here. There are also some useful files listed at the bottom of this page and directions on where to obtain them. I've not done any posts on accessing and processing the data (I code in Visual Basic), but if anyone wants me to cover that just ask and I will do so. While I'm on the subject of requests, I have been using gridded Drift Age Model data, although I didn't find it interesting enough to really get into, if anyone wants to know how to access that data just ask.

First the grid box statistics and effective thickness.
Latitude:      84.12
Longitude:     15.36
Area:          580.7655
Region:        10
EffThick:      1.992985 
The grid box is at the stated latitude and longitude near the pole, it's in region 10 (Central Arctic), and has an area of 580.77km^2. The effective thickness of the grid box is 1.993m, but that isn't the whole story regards thickness.

Now for the gice data, this is presented as the centre thickness of the thickness category (in metres), and then the gice thickness distribution for this grid cell. Those latter numbers are printed as they come out of the file, the printed version of single precision floating point numbers (IEE 32 bit, 4 bytes per number) .
 0             0.1388317
 0.26          2.307221E-02
 0.71          9.981988E-02
 1.46          0.3304134
 2.61          0.2930775
 4.23          8.014728E-02
 6.39          1.488431E-02
 9.1           8.990611E-03
 12.39         6.379645E-03
 16.24         4.273801E-03
 20.62         9.55024E-05
 25.49         1.397002E-05 
These gice numbers are graphed in the following column plot. Note that they are a distribution normalised to 1, in other words the sum of all twelve gice numbers is 1.

Looking at the above plot it can be seen that even in the centre of the pack in December there is a small but significant presence at 0m, in other words open water in leads between floes. Amounts in the thinner categories are small, and the gice distribution peaks at around 1.46m and 2.61m. the presence of non-zero gice numbers for ice 12.39m and 16.24m shows the presence of ridged ice, indeed there is a small presence of thicker ice within the grid box.

But how does this relate to effective thickness? That starts by multiplying each gice number with its respective thickness, which generates a new series of numbers (truncated here for clarity).

0 0.0000
0.26 0.0060
0.71 0.0709
1.46 0.4824
2.61 0.7649
4.23 0.3390
6.39 0.0951
9.1 0.0818
12.39 0.0790
16.24 0.0694
20.62 0.0020
25.49 0.0004

In this set however the sum of the numbers is not 1, it is 1.990931, that is the grid box effective thickness. In practice it is a fraction of a percent too low. I cannot identify why it isn't a closer match, my calculation of volume using gridded heff files has an average error of 4.74ppm and a standard deviation of 71.14ppm. Volume calculated for December using gice has an error of -0.17% average, with a standard deviation of 0.034%. I've emailed Dr Zhang, but haven't been able to get to the bottom of the matter. So I will continue to use heff to calculate regional area, while using gice to get a look at the sub grid thickness distribution and the presence of thicker and thinner categories of ice. See comments for the reason for this.

So staying with this example grid box, what does the thickness distribution look like?

While the effective thickness of 1.993m might imply a homogeneous layer of about 2m thick, the distribution above shows that within PIOMAS the picture is not so simple. By definition the open water seen in the gice distribution plays no role, any gice number multiplied by 0 thickness is 0. Now the peak is at 2.61m and quite clearly so. But all this is getting rather confusing for me. The sum of those numbers is about 2m, so each of those figure is a fraction of 2m, so we're saying that the 12.39m thickness band is 7.9cm of the total 2m... frankly I don't really know how to get my head around that. I feel like I'm following a white rabbit who's in a hurry.

There is a better way to get a handle on thickness distribution. Going back to how the effective thickness is calculated from the gice numbers we can use area to introduce volume for each thickness band, which is just the thickness of the gice band multiplied by its respective gice number multiplied by the grid cell area.

For this representation the sum of the volumes for the different bands roughly equals the effective thickness multiplied by the grid box area (1.1575km^3), and exactly equals the equivalent calculation using gice (1.1563km^3). But it gives a simple mind like mine more of a grasp of what is going on. For example the volume for the 6.39 to 16.24m bands is about the same as for the 4.23m band, that's a significant amount of very thick ice. Much easier to grasp than 12.39m thickness band making up 7.9cm of the total 2m effective thickness.

I will be incorporating gice into my regular posts, but will also have a look at certain details in the weeks to come. For the present, I've graphed gice for the Arctic Ocean in September showing the impacts of the losses of 2007, 2010 and 2012.

In that graph the colour bands follow the same order as the key on the right hand side of the plot.


Chris Reynolds said...

I've been discussing with Dr Zhang the difference between volume calculated from heff and volume calculated from monthly gice. Also the problem of calculating monthly gice from daily.

Dr Zhang has just realised that as daily gice is calculated at the end of the day, while heff (effective thickness) uses all gice during the day, using all 75 timesteps in the day for each day of the month. Whereas monthly mean gice uses all 75 time steps for each day of the month.

The further problem I have been having is that my calculation of monthly average gice isn't matching up with the monthly gice data. For this reason I will be using monthly gice for long term work. But for comparing state in 2015 (for which monthly gice is not yet available) I will have to use the daily data. So in April state will be calculated using 30 April data.

Brent Ragsdale said...


I enjoy reading your blog. Thanks for posts like this one that attempt to explain the analysis in plain English. A while back you wrote about a mono-modal spike. That was quite thought provoking for me. I wonder if the implications of your main graph in that post don't provide a possible explanation for why the arctic ice loss might be slowing down such that the Gompertz curve seems like the best fit. If you extrapolate that graph into the future I wonder if it will play out something like this:
Anyway, I think you are on to something by breaking down grid thickness. Cheers! Brent Ragsdale-Kansas City

Chris Reynolds said...

Thanks for that graphic Brent.

It's a very impressive visualisation and I think captures what we're going to see (on average) in the decades to come.

You've illustrated the decline in the peak of thickness as ice thins with warmer winters, and the resulting decline in peak volume. The 'long tail' of thickness towards the higher thickness is also to be expected to persist due to ridging of ice and survival of ice off the Canadian Arctic Archipelago.

It will be interesting to see later this century if volume will decline as your illustration suggests or if there will be a rapid collapse. Some of the research suggests that while there is no 'bifurcation' in the transition to a seasonally sea ice free state, there may be such a non-linear collapse in the transition to perennially ice free. e.g. Eisenman "Factors controlling the bifurcation structure of sea ice retreat."

It's not me that is onto something with breaking down grid thickness. Dr Zhang is aware of this blog and the use Wipneus has made of the gridded data, and Dr Zhang made this sub grid thickness data available some time ago. I have actually been rather slow in getting round to looking at it. I suspect Dr Zhang might have been prompted by my thickness distribution treatment of effective thickness - which is useful but doesn't tell the whole story.