Lake ice thickness

OOPS: I posted this late at night and as I laid down to bed afterwards I realized that I was using the wrong \sigma. I was using the one for water, not ice. Ice is larger by a factor of 4, so that explains why my estimate was so low. I’ll redo the calcs hopefully later today (basically I just multiply everything by 4).

This year my partner gave me a manual ice auger for Christmas. I love fishing, and I thought it might be fun to ice fish as well. I do enjoy it, but it’s definitely different.

Here in Minnesota, we get reports every year of people falling through the ice. Sometimes just themselves, sometimes with their snowmobiles or their trucks. When I go, I usually park in the lots and walk, but I will admit I’m tempted to drive on out there.

So I’ve been thinking about how to determine the thickness of ice. I thought I’d take a simple approach and see how far I could get. It turns out, you can make it quite a ways, and along the way I discovered some interesting issues about numeric calculations.

Underlying physics

Here’s the way I approached it: We need to determine how long it’ll take to freeze the next layer of water that’s just under the existing ice. I often tell my students that to solve growth problems like these, you should jump to the middle, figure out the governing physics for the next time step, and then integrate with the appropriate boundaries: in this case, thickness=0 at the beginning.

So, let’s assume we have some thickness of ice, x(t). The next layer of ice, of thickness dx, will freeze once you’ve gotten the appropriate amount of energy sucked out of it, through the ice, into the cold air right above the ice.

dQ=m L

Where dQ is the amount of energy you’re extracting, m is the mass of the next layer of ice, and L is the heat of fusion (334,000 J/kg) of water. The mass of that narrow slice is given by:

m=\rho A dx

where A is the area of the lake (don’t worry, as I always tell my students, if something you don’t know creeps in, watch for it to cancel).

Ok, how fast can we suck that heat through the existing ice? For that we need the formula for conduction. However, I always like to think of it as Ohm’s law (V=IR) where V is replaced by the temperature difference, I is replaced by the heat flow (dQ/dt), and R is replaced by the heat resistance:

R=\frac{\text{thickness}}{\sigma A}

where \sigma is the heat conductivity or 0.5985 W/m-K.

So, let’s put it all together:

\Delta T=\frac{x}{\sigma A}\frac{dQ}{dt}

\Delta T=\frac{x}{\sigma A}\rho A L\frac{dx}{dt}

\Delta T=\frac{\rho L}{\sigma}x \frac{dx}{dt}

Note how the Areas cancelled (just like I told you they would). For \Delta T we assume that the temperature below the ice is right at freezing and the temperature above the ice we get from weather stations. If \Delta T is positive, the ice thickness grows, so we want T_below – T_above.

Ok, now we’ve got a differential equation. Note that the temperatures are functions of time, so we have to be careful. Here’s a couple of ways to go about it: 1) Stick that differential equation into a differential equation numerical solver (with the temperatures interpolated) and plot the result, and 2) analytically integrate the equation to get the result from a direct integration of the temperatures. I’ll give away a portion of the punch line right now: First, (1) and (2) give the same results. Second, (1) is way faster than (2) in Mathematica, though, for the moment, I’m not 100% sure why.

Differential equation solver

My longtime readers will surely know that NDSolve is coming off the bench, once again. I grab the temperature data from a nearby weather station, interpolate it so that Mathematica can calculate the air temperature at any time I like in the month of December, and then I put the above differential equation into NDSolve and see what I get.

CAVEAT: You have to start all of these calculations at the right time of the month. We didn’t get cold and stay cold until about 10 days in to December this year. I look at the temperature data and pick the day when we went freezing and stayed freezing for a while. I also need to start with a tiny amount of ice, or else there’s a divide by zero that Mathematica always wants to do.

By the way, you can get the temperature data right in Mathematica with a command like WeatherData["KSTP", "Temperature", {2012, 12}] where “KSTP” is the name of a local weather station. To figure that out, I actually tried “= temperature in st. paul in december 2012″ to activate the connection to Wolfram Alpha. Mathematica sends your request to WA and you get lots of stuff back (basically the same stuff WA would give you). You can grab the data from there, but the WeatherData command puts it in a better form. Here’s a typical result from WA:

temperature in St. Paul, MN in December 2012 from Wolfram Alpha

temperature in St. Paul, MN in December 2012 from Wolfram Alpha

Ok, I grab the data, interpolate it, and put in my NDSolve command. Here’s the result:

Lake ice thickness (in meters) versus time in December (seconds)

Lake ice thickness (in meters) versus time in December (seconds)

What’s weird is that I went ice fishing at the end of the month, and I drilled for at least twice that thickness. Not sure what to make of that.

Analytical integration

It turns out the differential equation in question is not that hard to solve. First, bring the dt over to the left hand side, and then integrate both sides:

\int \Delta T(t)\,dt = \frac{\rho L}{\sigma}\frac{1}{2}x^2

OOPS: I originally had the constants outside of the square root. It’s fixed now.

x(t)=\sqrt{\frac{2\sigma}{\rho L}\int \Delta T(t)\,dt}

where, of course, you have to be careful of the limits of the integral. For x, we assume the ice starts at zero thickness. For time, we start on the same day mentioned above, and we integrate to whatever time we’d like.

What’s really cool about this approach, is you can see that if the integral is zero, x is zero. In other words, if you have just as much area above freezing as below, the ice will be gone. You can also directly see how \sigma, \rho and L affect the result. Higher conductivity? more ice (because it can suck the heat out faster). Higher density? less ice (more energy to suck out). Higher heat of fusion? less ice (more energy to suck out).

What’s weird about this approach, is to do the integration of the time series data you get from Wolfram Alpha, I tried to simply integrate the interpolation function I made of the temperatures. This took something like 30 times as long as the NDSolve command. I’m not really sure why. When I wrote my own integration function for the time series data, it was lickety-split, just like NDSolve. Weird.

The upshot, though, is that I get the same result.

Which do you like better?

About these ads

About Andy "SuperFly" Rundquist

Associate professor of physics at Hamline.
This entry was posted in mathematica, physics. Bookmark the permalink.

4 Responses to Lake ice thickness

  1. ntmoore says:

    Wait, how is the energy getting out of the ice? I wonder if the results would better match the ice depth if you added some sort of convective (or radiative) heat transfer that removes heat from the surface layer.

    • Andy "SuperFly" Rundquist says:

      I guess I was assuming that once the heat gets to the top, it leaves much faster than the transfer through the ice. With my new correction (see new first paragraph), I’m actually right around (or maybe a little higher than) what I experienced on the lake.

  2. Your factor of 4. Ice has four modes for phonons vibrations to travel torsion, pressure, two shear modes. Water only has a pressure mode. 1+1+2 modes / 1 mode ~ 4x (I overlooked the torsional mode in my tweet)

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s