Thermomechanisches Verhalten polythermer Eisschilde - Theorie, Analytik, Numerik

Thermomechanisches Verhalten polythermer Eisschilde – Theorie, Analytik, Numerik – Vom Fachbereich Mechanik der Technischen Hochschule Darmstadt zur Erlangung des akademischen Grades eines Doktors der Naturwissenschaften genehmigte Dissertation von Dipl.-Phys. Ralf Greve aus Siegburg Referent: Prof. K. Hutter, Ph. D./Cornell Univ. Korreferent: Prof. Dr. K. Herterich Tag der Einreichung: 11.5.1995 Tag der m¨ undlichen Pr¨ ufung: 30.8.1995 Darmstadt 1995 D 17 Here is a door ajar through which one may escape a little way and for a short time out of our little world, from the noise and chaos of civilization into the silence and harmony of the cosmos, and for a moment be part of it. – Richard E. Byrd Danksagung Diese Dissertation wurde durch Herrn Prof. Kolumban Hutter, Ph. D., angeregt und erm¨ oglicht. Ich danke ihm f¨ ur die intensive Unterst¨ utzung und F¨ orderung, welche mir von seiner Seite w¨ ahrend der Anfertigung der Arbeit zuteil wurden, sowie f¨ ur die sorgf¨ altige Durchsicht des Manuskripts. Des weiteren m¨ ochte ich Herrn Prof. Dr. Klaus Herterich f¨ ur die ¨ Ubernahme des Korreferates danken. Ich danke meinen Kollegen und Ex-Kollegen der Arbeitsgruppe III f¨ ur die Zusammenarbeit und das stets gute Arbeitsklima. Insbesondere sei den Herren Dr. Stefan Diebels, Dipl.-Phys. Georg Bauer und Magnus Weis f¨ ur ihre zeitaufwendige Betreuung der institutseigenen Workstations, Herrn Dr. Reinhard Calov f¨ ur viele Diskussionen und seine Hilfe bei den Gr¨ onland-Simulationen, sowie Frau Dipl.-Oz. Imke Hansen f¨ ur gr¨ undliches Korrekturlesen gedankt. Herrn Dipl.-Ing. J¨ org Schneider von der Firma Technosystem danke ich f¨ ur die Hilfestellung bei der Betreuung der Macintosh-Computer am Institut sowie bei der Anschaffung und Unterhaltung eines privaten Ger¨ ates. Im Sommer 1994 hatte ich die M¨ oglichkeit, drei Monate am Department of Geophysical Sciences der University of Chicago zu arbeiten. Hierf¨ ur sei Herrn Prof. Douglas R. MacAyeal, Ph. D., gedankt. Der Studienstiftung des deutschen Volkes, und hier speziell meinem Darmst¨ adter Vertrauensdozenten, Herrn Prof. Dr. Achim Richter, und meinem Ansprechpartner in Bonn, Herrn Dr. Ulf Lange, danke ich f¨ ur die zweieinhalbj¨ ahrige F¨ orderung dieser Dissertation. Weiterhin danke ich Frau Dr. Anne Letr´eguilly f¨ ur die Bereitstellung der Topographiedaten f¨ ur Gr¨ onland, sowie Herrn Dr. Sigfus Johnsen f¨ ur die ¨ Uberlassung der δ18O-Daten des GRIP-Eisbohrkerns. Schließlich gilt mein Dank meiner lieben Freundin Frau Anke Kurzke und meinen Eltern, Herrn Dieter und Frau Ingrid Greve. Ohne die von dieser Seite erfahrene moralische und finanzielle Unterst¨ utzung w¨ are die Vollendung dieser Arbeit nicht m¨ oglich gewesen. Inhaltsverzeichnis Zusammenfassung 5 Abstract 7 Notation 9 1 Einleitung 15 1.1 Klimageschichte der Erde . . . . . . . . . . . . . . . . . . . . . . . . . 15 1.2 Modellierung von Inlandeisschilden . . . . . . . . . . . . . . . . . . . 17 2 Kontinuumsmechanische Grundlagen 23 2.1 Allgemeine Bilanzaussagen . . . . . . . . . . . . . . . . . . . . . . . . 23 2.2 Massenbilanz ............................... 24 2.3 Impulsbilanz................................ 24 2.4 Drehimpulsbilanz............................. 25 2.5 Energiebilanz ............................... 25 2.6 Entropiebilanz............................... 25 2.7 Materialgesetze f¨ urviskoseFluide.................... 26 3 Das polytherme Eismodell 29 3.1 Feldgleichungen.............................. 30 3.1.1 KalterBereich........................... 30 3.1.2 Temperierter Bereich . . . . . . . . . . . . . . . . . . . . . . . 31 3.1.3 Lithosph¨ are ............................ 34 3.2 Randund ¨ Ubergangsbedingungen . . . . . . . . . . . . . . . . . . . 35 3.2.1 Randbedingungen an der freien Oberfl¨ ache........... 35 3.2.2 ¨ Ubergangsbedingungen an der kalten Eisbasis . . . . . . . . . 37 3.2.3 ¨ Ubergangsbedingungen an der temperierten Eisbasis . . . . . 38 3.2.4 Randbedingungen an der Lithosph¨ arenunterseite . . . . . . . . 41 3.2.5 ¨ Ubergangsbedingungen an der CTS . . . . . . . . . . . . . . . 41 4 Planparallele, geneigte Eisplatte 47 4.1 Anwendung des Modells . . . . . . . . . . . . . . . . . . . . . . . . . 47 4.2 Integration der Slab-Gleichungen . . . . . . . . . . . . . . . . . . . . 50 4.3 Ergebnisse................................. 55 1 5 Skalierung und Flacheisannahme 59 5.1 Einf¨ uhrungderSkalierung ........................ 59 5.2 Skalierung und SIA f¨ ur die Modellgleichungen . . . . . . . . . . . . . 62 5.2.1 KalterBereich........................... 62 5.2.2 Temperierter Bereich . . . . . . . . . . . . . . . . . . . . . . . 64 5.2.3 Lithosph¨ are ............................ 65 5.2.4 Randbedingungen an der freien Oberfl¨ ache........... 66 5.2.5 ¨ Ubergangsbedingungen an der kalten Eisbasis . . . . . . . . . 67 5.2.6 ¨ Ubergangsbedingungen an der temperierten Eisbasis . . . . . 70 5.2.7 Randbedingungen an der Lithosph¨ arenunterseite . . . . . . . . 71 5.2.8 ¨ Ubergangsbedingungen an der CTS . . . . . . . . . . . . . . . 71 5.3 Teilintegration der polythermen SIA-Gleichungen . . . . . . . . . . . 73 5.3.1 Berechnung der Spannungen . . . . . . . . . . . . . . . . . . . 74 5.3.2 Berechnung der Geschwindigkeit . . . . . . . . . . . . . . . . . 74 5.3.3 Evolution der freien Oberfl¨ ache ................. 76 5.3.4 Evolution der CTS . . . . . . . . . . . . . . . . . . . . . . . . 76 5.3.5 Evolution der Eisbasis bzw. Lithosph¨ arenoberseite . . . . . . . 77 5.3.6 Temperatur und Wassergehalt . . . . . . . . . . . . . . . . . . 77 5.4 Zusammenstellung in dimensionsbehafteter Form . . . . . . . . . . . 80 5.5 Spezifizierung physikalischer Gr¨ oßen................... 83 6 Numerische L¨ osung der polythermen SIA-Gleichungen 87 6.1 σ-Transformation und Transformationsregeln . . . . . . . . . . . . . . 87 6.2 Transformierte Modellgleichungen . . . . . . . . . . . . . . . . . . . . 88 6.3 Das numerische Gitter . . . . . . . . . . . . . . . . . . . . . . . . . . 92 6.4 Diskretisierung der Modellgleichungen . . . . . . . . . . . . . . . . . . 94 6.5 Positionierung der CTS . . . . . . . . . . . . . . . . . . . . . . . . . . 108 6.6 Das Programmpaket SICOPOLIS . . . . . . . . . . . . . . . . . . . . 116 7 Simulationen f¨ ur das EISMINT-Eisschild 119 7.1 Das EISMINT-Eisschild . . . . . . . . . . . . . . . . . . . . . . . . . 119 7.2 Diskussion der Simulationen . . . . . . . . . . . . . . . . . . . . . . . 121 7.3 Abbildungen zu den Simulationen f¨ ur das EISMINT-Eisschild . . . . 127 8 Simulationen f¨ ur das gr¨ onl¨ andische Eisschild 153 8.1 Stereographische Projektion . . . . . . . . . . . . . . . . . . . . . . . 153 8.2 Die heutige Topographie Gr¨ onlands................... 154 8.3 Standard-Randbedingungen f¨ ur heutiges Klima . . . . . . . . . . . . 155 2 8.4 Standard-Randbedingungen f¨ ur andere Klimaszenarien . . . . . . . . 157 8.5 Diskussion der Simulationen . . . . . . . . . . . . . . . . . . . . . . . 158 8.5.1 Steady-State-Simulationen . . . . . . . . . . . . . . . . . . . . 158 8.5.2 Transiente Simulationen mit Milankovi´c-Sinus-Antrieb . . . . 168 8.5.3 Transiente Simulationen mit GRIP-Daten-Antrieb . . . . . . . 169 8.5.4 Transiente Simulationen zur Auswirkung eines anthropogen verst¨ arkten Treibhauseffektes . . . . . . . . . . . . . . . . . . . . 175 8.6 Abbildungen zu den Simulationen f¨ ur das gr¨ onl¨ andische Eisschild . . 177 9 Ausblick 219 Literaturverzeichnis 221 3 gSchwerebeschleunugung h z-Koordinate der freien Eisoberfl¨ ache hmax maximales h des gesamten Eisschilds HEisdicke Hmax maximale Eisdicke H HcDicke der kalten Eisschicht HtDicke der temperierten Eisschicht Ht,max maximale Dicke Htder temperierten Eisschicht Hoffset terstes probeweise angenommenes Htbei neu entstehender temperierter Schicht Ht,smooth numerisch gegl¨ attetes Htvor der Korrektur zur Volumenerhaltung HrDicke der Lithosph¨ are i,j,kc,t,r Laufindices f¨ ur die diskretisierten x-, ybzw. ζc,t,r-Richtungen imax,jmax Maximalwerte f¨ ur i,j jdiffusive Wasserstromdichte jtot totale Wasserstromdichte kc,t,r,max Maximalwerte f¨ ur kc,t,r KStreckfaktor f¨ ur die stereographische Projektion Llatente Schmelzw¨ arme f¨ ur Eis Mi) Wasserproduktionsrate im temperierten Eis ii) Ablationsrate (Schmelzrate an der Eisoberfl¨ ache) M?Bildungsrate f¨ ur ¨ Uberlagerungseis ˙mw bWasser-Massenstrom in den Boden hinein n, ˜nLaufindices f¨ ur die feinere bzw. grobere Diskretisierung der Zeit nmax, ˜nmax Maximalwerte f¨ ur n, ˜n pDruck Pmax S¨ attigungsgrad von Wasser in Schnee Q⊥ geoth geothermer W¨ armefluß qW¨ armefluß (nur f¨ uhlbare W¨ arme) qges W¨ armefluß (f¨ uhlbare und latente W¨ arme) q,qihorizontaler Massenfluss, bzw. dessen i-Komponente Runiverselle Gaskonstante ReErdradius SAkkumulationsrate (Schneefallrate) tZeit tacc Akkumulations-Zeitpunkt f¨ ur ein Eispartikel tend Abbruch-Zeit f¨ ur numerische Simulationen 10 tMil Periode der Milankovi´c-Zyklen TTemperatur TMdruckkorrigierte Schmelztemperatur T0Schmelztemperatur bei verschwindendem Druck T0homologe Temperatur (T−TM) T0 div T0in der S¨ aule (xdiv, ydiv) T0 half T0in der S¨ aule (xhalf, ydiv) Tair Lufttemperatur ¨ uber dem Eis ˆ Tair Amplitude der j¨ ahrlichen Variation von Tair Tma Jahresmittel der Lufttemperatur Tair TbTemperatur an der Eisbasis T0 bhomologe Temperatur an der Eisbasis T0 b,div T0 bin der S¨ aule (xdiv, ydiv) TsEisoberfl¨ achentemperatur (10 m-Firntemperatur) TSpannungstensor TRSpannungsdeviator (Reibungsspannungen) v,viGeschwindigkeit (baryzentrisch), bzw. deren i-Komponente vb, (vb)iEisgeschwindigkeit (baryzentrisch) an der Basis, bzw. deren i-Komponente vs, (vs)iEisgeschwindigkeit (baryzentrisch) an der Oberfl¨ ache, bzw. deren i-Komponente vsl, (vsl)iBasale Gleitgeschwindigkeit (Differenz zwischen Eisgeschwindigkeit und Lithosph¨ arengeschwindigkeit), bzw. deren i-Komponente vi,vi,j Geschwindigkeit des Eises in der Mischung Eis plus Wasser, bzw. deren j-Komponente vw,vw,j Geschwindigkeit des Wassers in der Mischung Eis plus Wasser, bzw. deren j-Komponente vx,half vxin der S¨ aule (xhalf, ydiv) vz,div vzin der S¨ aule (xdiv, ydiv) vz,half vzin der S¨ aule (xhalf, ydiv) Vges gesamtes Eisvolumen Vtemp Volumen des temperierten Eises Vsmooth temp Volumen des temperierten Eises nach der numerischen Gl¨ attung von Ht w,wiGeschwindigkeit einer singul¨ aren Fl¨ ache, bzw. deren i-Komponente x,yhorizontale kartesische Koordinaten xdiv,ydiv Position (x, y) des Scheitelpunktes xhalf x-Koordinate auf halbem Weg zwischen Scheitelpunkt und Eisrand 11 x2x-Koordinate des vom Scheitelpunkt aus gesehen zweiten Gitterpunktes mit temperierter Eisschicht l¨ angs ydiv zvertikale kartesische Koordinate (H¨ ohe ¨ uber NN) zmz-Koordinate der CTS zshift mBetrag der probeweisen CTS-Verschiebungen bei deren Positionierung αVerh¨ altnis potentielle Energie zu innere Energie in kaltem Eis αtVerh¨ altnis potentielle Energie zu innere Energie in temperiertem Eis βClausius-Clapeyron-Gradient β1,2Koeffizienten f¨ ur Schmelz-Parameterisierung ∆x, ∆yOrtsaufl¨ osung f¨ ur die xund y-Richtung ∆ζc,t,r Ortsaufl¨ osung f¨ ur die ζc,t,r-Richtung ∆t,f ∆tfeinere bzw. grobere zeitliche Aufl¨ osung ∆Tma(t) Abweichung des Jahresmittels der Lufttemperatur ¨ uber dem Eis Tma von den heutigen Werten ∆Ts(t) Abweichung der Eisoberfl¨ achentemperatur Tsvon den heutigen Werten εi) innere Energie ii) Aspektverh¨ altnis κW¨ armeleitf¨ ahigkeit f¨ ur Eis κrW¨ armeleitf¨ ahigkeit f¨ ur die Lithosphere λgeographische L¨ ange λ0L¨ angenkreis, der unter der stereographischen Projektion auf die y-Achse abgebildet wird µfwc Firnerw¨ armungs-Koeffizient νWasser-Diffusivit¨ at im temperierten Eis ρwahre Dichte f¨ ur Eis bzw. die Mischung Eis plus Wasser ρiPartialdichte f¨ ur Eis in der Mischung Eis plus Wasser ρwPartialdichte f¨ ur Wasser in der Mischung Eis plus Wasser ρawahre Dichte der Asthenosph¨ are ρrwahre Dichte der Lithosph¨ are σeffektive Scherspannung σij ij-Komponente des Spannungstensors σR ij ij-Komponente des Spannungsdeviators τVVerz¨ ogerungszeit der isostatischen Einstellung der Asthenosph¨ are τc,t,r Zeit im σ-Koordinatensystem f¨ ur kaltes Eis, temperiertes Eis bzw. die Lithosph¨ are ξc,t,r,ηc,t,r horizontale Koordinaten des σ-Koordinatensystems f¨ ur kaltes Eis, temperiertes Eis bzw. die Lithosph¨ are 12 ζc,t,r vertikale Koordinate des σ-Koordinatensystems f¨ ur kaltes Eis, temperiertes Eis bzw. die Lithosph¨ are φgeographische Breite φ0Breitenkreis der Bildebene der stereographischen Projektion ωWassergehalt des temperierten Eises (Massenanteil) ωmax Schwellwert im Wasserdrainage-Gesetz ω2ωin der S¨ aule (x2, ydiv) [A] typischer Rate-Faktor [c] typische spezifische W¨ arme f¨ ur Eis [cr] typische spezifische W¨ arme f¨ ur die Lithosph¨ are [C] typische Gleitfunktion f¨ ur kalte Basis [Ct] typische Gleitfunktion f¨ ur temperierte Basis [f] typische Kriechfunktion [H] typische Vertikaldimension [L] typische Horizontaldimension [Q⊥ geoth] typischer geothermer W¨ armefluß [VH] typische Vertikalgeschwindigkeit [VL] typische Horizontalgeschwindigkeit [∆T] typische Temperaturdifferenz [κ] typische W¨ armeleitf¨ ahigkeit f¨ ur Eis [κr] typische W¨ armeleitf¨ ahigkeit f¨ ur die Lithosph¨ are [ω] typischer Wassergehalt BClausius-Clapeyron-Zahl DW¨ armediffusionszahl DtWasserdiffusionszahl DrW¨ armediffusionszahl der Lithosph¨ are E(ζc) Streckfunktion f¨ ur die σ-Transformation des Kalteisbereichs FGleitzahl f¨ ur kalte Basis FtGleitzahl f¨ ur temperierte Basis KFluidit¨ atszahl Nrgeotherme W¨ armezahl in der Lithosph¨ are Pw bbasale Schmelzrate Pw mOberfl¨ achenproduktionsrate von Wasser an der CTS TrZeitverz¨ ogerungszahl f¨ ur die Asthenosph¨ are FFroude-Zahl 13 14 1 Einleitung 1.1 Klimageschichte der Erde Seit sich die Erde vor etwa 4.6 Milliarden Jahren aus einem kleinen Teil des kosmischen Urnebels, der unser Sonnensystem geworden ist, gebildet hat, unterlag ihr Klima einem best¨ andigen Wandel. Die ersten hundert Millionen Jahre waren als Folge der in W¨ arme umgewandelten Gravitationsenergie und der W¨ armeproduktion radioaktiver Isotope der neu entstandenen Erde sehr heiß mit Oberfl¨ achentemperaturen ¨ uber 100◦C, bis diese Grenze vor ungef¨ ahr vier Milliarden Jahren unterschritten wurde. Die Abk¨ uhlung der Erdoberfl¨ ache zog sich anschließend noch weitere zwei Milliarden Jahre hin, bis ca. zwei Milliarden Jahre vor heute ein relatives Gleichgewicht erreicht war; zu dieser Zeit hatte sich auch bereits eine Stickstoff-Sauerstoff-Atmosph¨ are, ¨ ahnlich der heutigen, gebildet. Blickt man auf dieser Jahrmilliarden umfassenden Zeitskala in die Zukunft, so wird sich die Erdoberfl¨ ache aufgrund der zunehmenden Energieproduktion der Sonne langsam erw¨ armen und schließlich in ferner Zukunft wieder den Glutofen darstellen, der sie in der Fr¨ uhzeit der Erdgeschichte einmal war. Dies beschreibt jedoch nur den langfristigen, die gesamte Lebensdauer der Erde umfassenden Trend der Klimageschichte. Diesem Trend ¨ uberlagern sich ¨ außerst vielf¨ altige Variationen, die sich auf Zeitskalen von 100 Millionen Jahren bis hin zu Jahrzehnten abspielen. (Nat¨ urlich existieren auch noch kurzfristigere Schwankungen, diese werden allerdings nicht mehr als Klima, sondern als Wetter bezeichnet.) Bleibt man zun¨ achst bei den großen Zeitskalen von 10-100 Millionen Jahren, so zeigt sich, daß es in der Klimageschichte der Erde Wechsel zwischen nicht-eisbildenden Warmklimazeitaltern einerseits und Eiszeitaltern andererseits gegeben hat. Warmklimazeitalter zeichnen sich dadurch aus, daß die Erdoberfl¨ ache v¨ ollig eisfrei war, w¨ ahrend in Eiszeitaltern Eisbedeckungen auftraten. Diese Phasen sind sehr ungleich verteilt; auch bei Beschr¨ ankung auf die letzten 2.3 Milliarden Jahre der irdischen Klimageschichte seit dem Auftreten des ersten Eiszeitalters (“Archaisches Eiszeitalter”) ¨ uberwiegen die warmen Phasen die sieben identifizierten Eiszeitalter deutlich an Dauer. Die haupts¨ achliche Ursache f¨ ur das Auftreten von Eiszeitaltern ist wahrscheinlich das Vorhandensein von ausgedehnten Landmassen in polaren Regionen aufgrund der Kontinentaldrift, denn nur dann ist die Bildung großer Eismassen aus akkumuliertem Schneefall in polaren Regionen ¨ uberhaupt m¨ oglich. In diesem Fall tritt eine starke positive R¨ uckkopplung auf, die ¨ uber die hohe Albedo (R¨ uckstrahlverm¨ ogen des einfallenden Sonnenlichts) der vereisten Fl¨ achen zu verringerter Energieaufnahme der Erdoberfl¨ ache und damit zu Abk¨ uhlung f¨ uhrt. Sind die polnahen Gebiete dagegen haupts¨ achlich von Ozeanen bedeckt, wird eventuell fallender Schnee schnell geschmol15 zen, so daß keine großfl¨ achige Vereisung m¨ oglich ist. Nat¨ urlich setzt diese Erkl¨ arung voraus, daß das Klima hinreichend kalt ist, um Schneefall in polaren Regionen zuzulassen; daher gab es in den ersten 2.3 Milliarden Jahren nach Entstehung der Erde wahrscheinlich keine Eiszeitalter. W¨ ahrend des gesamten Erdmittelalters (Mesozoikum, 225-65 Millionen Jahre vor heute, bestehend aus Trias, Jura und Kreide), nachdem die Erde bereits sechs identifizierte Eiszeitalter erlebt hatte, herrschte eine Warmklimaphase vor, in der viele nette kleine Tierchen, genannt Dinosaurier, auf der Erde umherliefen. Nach dem Beginn des Terti¨ ars der Erdneuzeit (Neozoikum, seit 65 Millionen Jahren vor heute) setzte eine stufenweise Abk¨ uhlung ein, einhergehend mit der Bewegung des antarktischen Kontinents in die s¨ udpolare Region und dessen Vereisung in der zweiten H¨ alfte des Terti¨ ars. Der eigentliche Beginn dieses sich herauskristallisierenden Eiszeitalters wird jedoch ¨ ublicherweise erst auf zwei bis drei Millionen Jahre vor heute, den ¨ Ubergang vom Terti¨ ar zum Quart¨ ar, datiert, zusammenfallend mit einer weiteren sprunghaften Abk¨ uhlung. Dieses bis heute und wahrscheinlich noch viele Millionen Jahre in die Zukunft (falls nicht durch menschliche Eingriffe verhindert) andauernde Eiszeitalter wird daher als quart¨ ares Eiszeitalter bezeichnet. Auf einer Zeitskala von 104-105Jahren ist das quart¨ are Eiszeitalter gepr¨ agt vom Wechselspiel von mindestens 20 Kaltzeiten (“Eiszeiten”) und dazwischenliegenden Warmzeiten. Diese Warmzeiten sind nicht zu verwechseln mit Warmklimazeitaltern; auch die Warmzeiten eines Eiszeitalters weisen ausgedehnte Eisbedeckungen auf, wie die heutige, seit etwa 11000 Jahren herrschende Neo-Warmzeit zeigt. Davor lagen die W¨ urm-Kaltzeit (ab ca. 70000 Jahren vor heute), die Eem-Warmzeit (ab ca. 125000 Jahren vor heute), die Riß-Kaltzeit (ab ca. 200000 Jahren vor heute) usw. Dieses Wechselspiel geschieht vor dem Hintergrund einer praktisch station¨ aren Verteilung der Landmassen und kann daher nicht mehr mit der Kontinentaldrift erkl¨ art werden. Die prim¨ are Ursache liegt statt dessen in periodischen Schwankungen von Parametern der Umlaufbahn der Erde um die Sonne (Orbitaltheorie; Milankovi´c [54], Hays, Imbrie & Shackleton [21], Berger [5]). Hierbei treten drei Perioden auf, n¨ amlich die Schwankung der Bahnexzentrizit¨ at mit einem Zyklus von 96000 Jahren, die Schwankung der Erdachsenneigung relativ zur Bahnebene (40000 Jahre) und die Pr¨ azession der Erdachse zusammen mit der Orientierung der Bahnellipse im Raum (21000 Jahre), welche die Verteilung der Sonneneinstrahlung auf die Erdoberfl¨ ache ver¨ andern und ¨ uber mannigfache R¨ uckkopplungen positiver und negativer Art diese Klimaver¨ anderungen bewirken. Auch innerhalb von Kaltund Warmzeiten existieren markante Klimaschwankungen. W¨ ahrend der letzten (W¨ urm-) Kaltzeit traten drei besonders kalte Epochen, 16 sogenannte Stadiale, vor ca. 60000, 40000 und 18000 Jahren auf (letztere ist bekannt als “Last Glacial Maximum”, oder kurz LGM), unterbrochen von w¨ armeren Interstadialen. Das Ausmaß von Klimaschwankungen in der Eem-Warmzeit davor ist aufgrund sich teilweise widersprechender Rekonstruktionen umstritten (Zahn [72]); die heutige (Neo-) Warmzeit ist gepr¨ agt von einem relativ gleichf¨ ormigen Klima mit Temperaturschwankungen im Bereich von maximal 1◦C seit etwa 10000 Jahren, was eventuell eine große Ausnahmesituation darstellt (Stauffer [67]), die eine Entwicklung der menschlichen Kultur (?) erst m¨ oglich gemacht hat. Jedoch gab es auch im Rahmen dieser Schwankungen ausgezeichnete Perioden, wie etwa das mittelalterliche Klimaoptimum von ca. 800-1200 n. Chr., die sich daran anschließende “kleine Eiszeit” (ca. 1400-1900 n. Chr., letzter Tiefpunkt um 1850), und schließlich das heutige, wieder etwas w¨ armere Klima, welches bis jetzt noch im Rahmen der letzten Jahrtausende liegt, sehr wahrscheinlich aber durch menschliche Eingriffe in das komplexe Klimasystem der Erde zumindest mitverursacht wurde und wird.1 1.2 Modellierung von Inlandeisschilden Die irdische Kryosph¨ are (der vereiste Teil des Klimasystems) setzt sich aus mehreren Komponenten zusammen. Unter einem Inlandeisschild, auch kurz als Inlandeis oder als Eisschild bezeichnet, versteht man eine ausgedehnte Eismasse, deren Basis auf festem Land aufliegt und die sich im Laufe von Jahrtausenden durch akkumulierten Schneefall gebildet hat. Inlandeisschilde machen bei weitem den gr¨ oßten Anteil des auf der Erde vorhandenen Eisvolumens aus; mit ihrer theoretischen und numerischen Modellierung besch¨ aftigt sich die vorliegende Arbeit. Alpine Gletscher entstehen analog, erstrecken sich jedoch nur ¨ uber wesentlich kleinr¨ aumigere Regionen in Gebirgen und tragen daher lediglich zu einem sehr geringen Teil zur Kryosph¨ are bei. Schelfeise sind aufschwimmende Eismassen, die vom meerw¨ artigen Massenfluß eines Inlandeises gespeist werden; sie existieren typischerweise in großen Buchten eines von Inlandeis bedeckten Kontinentalschelfs. Bei Meereis handelt es sich um gefrorenes Meerwasser, und Grundeis bedeutet schließlich gefrorenes Wasser im Boden, wie es in Permafrostgebieten auftritt. In Gletschern und Eisschilden vorkommendes Eis existiert (bei Vernachl¨ assigung zus¨ atzlicher Beimischungen von Salz und Sediment) in zwei grunds¨ atzlich verschiedenen Zust¨ anden. Kaltes Eis ist durch eine Temperatur unterhalb des druckkorrigierten Schmelzpunktes gekennzeichnet und kann als dichtebest¨ andiges, viskoses und w¨ armeleitendes Ein-Komponenten-Fluid beschrieben werden; es macht in den 1Die Ausf¨ uhrungen dieses Abschnitts folgen, falls nicht anders erw¨ ahnt, Sch¨ onwiese [66] und den dort angegebenen Prim¨ arquellen. 17 großen Eisschilden der Erde (Gr¨ onland, Antarktis) den weitaus gr¨ oßten Anteil aus. Bei temperiertem Eis dagegen befindet sich die Temperatur exakt auf dem druckkorrigierten Schmelzpunkt, was dazu f¨ uhrt, daß dieses zus¨ atzlich Wasser in geringer Menge enthalten kann. Demzufolge muß es im Gegensatz zu kaltem Eis als ZweiKomponenten-Fluid aufgefaßt werden. Bereiche aus temperiertem Eis k¨ onnen in Eisschilden in d¨ unnen, bodennahen Schichten existieren, die das Fließverhalten entscheidend ver¨ andern. Gletscher und Eisschilde, die sowohl kalte als auch temperierte Zonen enthalten, nennt man polytherm. In der Vergangenheit wurden bereits einige dreidimensionale Modelle zur numerischen Simulation des dynamischen Verhaltens von Inlandeisschilden entwickelt, was erst durch die hohen Rechenleistungen moderner Computer m¨ oglich gemacht wurde. Das erste Modell stammt von Mahaffy [53] und wurde angewendet auf die Barnes Ice Cap in der kanadischen Arktis. Es vernachl¨ assigte allerdings die starke Temperaturabh¨ angigkeit der Eisviskosit¨ at; so k¨ onnen die Modellgleichungen vertikal integriert werden, und das numerische Problem reduziert sich auf die zwei horizontalen Dimensionen. Das Modell von Jenssen [43] rechnete bereits mit temperaturabh¨ angiger Eisviskosit¨ at und wurde auf das gr¨ onl¨ andische Eisschild angewendet, allerdings mit sehr geringer r¨ aumlicher Aufl¨ osung aufgrund der damals noch beschr¨ ankten Rechnerkapazit¨ aten. Es folgten weitere Modelle, die auf unterschiedliche Probleme, wie Gr¨ onland, die Antarktis, Tibet, angewendet wurden (Oerlemans [58, 59], Budd & Smith [11], Budd, Jenssen & Smith [10], Herterich [23, 24], Kuhle, Herterich & Calov [45], Fastook & Chapman [16], Huybrechts & Oerlemans [40, 41], Huybrechts [34, 35, 36, 37, 38], Huybrechts, Letr´eguilly & Reeh [39], Letr´eguilly, Huybrechts & Reeh [47], Letr´eguilly, Reeh & Huybrechts [49], Calov [12], Calov & Hutter [13, 14], Fabr´e, Letr´eguilly, Ritz & Mangeney [15]). Herauszuheben sind die Antarktis-Simulationen von P. Huybrechts, welche ausf¨ uhrlich in [36] beschrieben sind; hierbei handelt es sich um das bis heute fortgeschrittenste Modell, welches das gekoppelte Inlandeis-Schelfeis-Lithosph¨ areProblem mit sehr hoher Aufl¨ osung rechnet. All diesen Modellen ist jedoch gemein, daß sie das m¨ ogliche Auftreten von temperierten Eisbereichen weitgehend vernachl¨ assigen. Modellrechnungen zur Bestimmung des Geschwindigkeitsund Temperaturfeldes von Eisschilden sowie der Evolution von Dicke und Ausdehnung bei gegebenem klimatischen Input werden mit den Prozeßgleichungen f¨ ur kaltes Eis durchgef¨ uhrt. Erh¨ alt man hierbei Temperaturen oberhalb des Schmelzpunktes, so wird f¨ ur die betreffenden numerischen Gitterpunkte die Temperatur k¨ unstlich auf diesen herabgesetzt. Diese Vorgehensweise stellt jedoch lediglich eine sehr grobe Ber¨ ucksichtigung der temperierten Bereiche dar, denn es wird nicht der Tatsache Rechnung getragen, daß mit den kalten und temperierten Bereichen zwei 18 verschiedene Phasen im Eisschild vorhanden sind, die durch eine Phasengrenzfl¨ ache voneinander getrennt sind, an welcher bestimmte Sprungbedingungen f¨ ur die verschiedenen physikalischen Gr¨ oßen erf¨ ullt sein m¨ ussen (M¨ uller [56], Hutter [29]). Ferner ist es nicht m¨ oglich, so den die Eisviskosit¨ at maßgeblich mitbestimmenden Wassergehalt in den temperierten Bereichen zu berechnen. Die besondere klimatologische Relevanz dieses Problems besteht darin, daß aufgrund der zu erwartenden, durch den verst¨ arkten Treibhauseffekt hervorgerufenen Erw¨ armung der Erdatmosph¨ are langfristig die große Gefahr des verst¨ arkten Abschmelzens der großen irdischen Eisschilde besteht, was einen Anstieg des Meeresspiegels um bis zu 70 m zur Folge h¨ atte. F¨ ur diesen Abschmelzprozeß sind basale temperierte Bereiche von vielleicht entscheidender Bedeutung. W¨ ahrend n¨ amlich kaltes Eis die Eigenschaft besitzt, an dem Felsgrund des Eisschildes zu haften, gehorcht temperiertes Eis dort aufgrund des Wassergehaltes einem viskosen Gleitgesetz; aus diesem Grunde k¨ onnen temperierte Bereiche das Abfließen der Eismassen (und somit deren Abbau) in die umgebenden Ozeane stark beschleunigen. Aus bisherigen Modellrechnungen mit Kalteismodellen und aus Radarmessungen l¨ aßt sich vermuten, daß große Bereiche der Antarktis basale temperierte Bereiche aufweisen (Huybrechts [36]); in Gr¨ onland existieren solche wahrscheinlich in der N¨ ahe des Eisrandes (Calov [12]). Diese Betrachtungen zeigen, daß es f¨ ur realistische Simulationen des Verhaltens großer Eisschilde sehr wichtig ist, die temperierten Bereiche ad¨ aquat und physikalisch korrekt zu behandeln. Ein weiterer Aspekt ist die Interpretation von Eisbohrkernen aus Gr¨ onland und der Antarktis. Hierbei geht man davon aus, daß das untersuchte Eis zu allen Zeiten kalt war, sich also nie Wasser bilden konnte, das durch das Eis diffundieren und so die Isotopenkonzentrationen und damit auch die Klimarekonstruktion verf¨ alschen konnte. Das Wissen um temperierte Bereiche und das Verhalten des Wassers darin k¨ onnte also eine ge¨ anderte Interpretation notwendig machen. Theoretische Beschreibungen zu polythermem Eis wurden bereits formuliert (Fowler & Larson [17], Hutter [27], Hutter [28], Blatter [6], Hutter [31]). Es werden zwei unterschiedliche Betrachtungsweisen verwendet: Diffusionsmodelle beschreiben temperiertes Eis mit zwei Massenbilanzen (eine f¨ ur die Mischung Eis plus Wasser, eine f¨ ur das Wasser selbst, dessen Anteil als gering angenommen wird), jedoch nur je einer Impulsund Energiebilanz f¨ ur die Mischung; der Transport des Wassers in der Eisumgebung wird dabei mit Hilfe eines Fickschen Diffusionsgesetzes dargestellt. Zwei-Phasenmodelle hingegen verwenden f¨ ur jede der beiden Komponenten sowohl eine Massenals auch eine Impulsbilanz, erg¨ anzt durch eine Wechselwirkungskraft zwischen den beiden Komponenten (Darcy-Transport). Erg¨ anzt werden diese Mo19 2.7 Materialgesetze f¨ ur viskose Fluide Das allgemeine Materialgesetz f¨ ur viskose Fluide (“Reiner-Rivlin-Fluide”) verkn¨ upft den Spannungstensor Tmit dem Verzerrungsgeschwindigkeits-Tensor Dund lautet T= (−p+ν0)1+ν1D+ν2D2(2.15) mit p=p(ρ, T), ν0,1,2=ν0,1,2(ρ, T, ID, IID, IIID),(2.16) ν0(D→0) = 0; hierbei ist pder Druck, ρdie Massendichte, Tdie Temperatur, ν0,1,2sind die Viskosit¨ atskoeffizienten, und ID,IID,IIIDbedeuten die drei Invarianten des Verzerrungsgeschwindigkeits-Tensors D(vgl. z. B. Hutter [30]). Im Falle eines dichtebest¨ andigen Fluids, wovon in der Eisdynamik ausgegangen wird, gilt gem¨ aß (2.6) div v= tr D= 0. Der Spannungstensor Twird in diesem Fall aufgespalten in einen isotropen Drucktensor und einen deviatorischen Tensor der Reibungsspannungen, T=−p1+TRmit tr TR= 0,(2.17) wobei der Druck jetzt ein freies Feld ist, und ein Materialgesetz nur noch f¨ ur TRzu formulieren ist. Dieses lautet also f¨ ur ein dichtebest¨ andiges Reiner-Rivlin-Fluid gem¨ aß (2.15) zun¨ achst TR=ν01+ν1D+ν2D2.(2.18) Durch Spurbildung erh¨ alt man ν0=−2 3ν2IID,(2.19) wobei 2 IID= tr D2verwendet wurde. (2.19) in (2.18) eingesetzt ergibt TR=ν1D+ν2D2−2 3IID1(2.20) mit ν1,2=ν1,2(T, IID, IIID) (2.21) als allgemeines Materialgesetz des dichtebest¨ andigen Reiner-Rivlin-Fluids. In der Eisdynamik wird ¨ ublicherweise die Vereinfachung getroffen, daß der nichtlineare Term proportional zu D2entf¨ allt, ν2also vernachl¨ assigbar ist. Weiterhin wird 26 der verbleibende Rest von (2.20) invertiert dargestellt. Da wegen Gl. (2.20), welche ja TRmit Dverkn¨ upft, die Invarianten von Ddurch die Invarianten von TRausgedr¨ uckt werden k¨ onnen, folgt D=1 ν1(T, IID, IIID)TR=: µ1(T, IITR,(IIITR)) TR.(2.22) Als weitere Vereinfachung wird, wie in (2.22) durch die Klammern bereits angedeutet, auch die Abh¨ angigkeit von der dritten Invarianten IIITRdes ReibungsspannungsTensors vernachl¨ assigt. Das verbleibende Materialgesetz ist jedoch im allgemeinen immer noch nichtlinear. 27 28 3 Das polytherme Eismodell Ein polythermes Eisschild besteht, wie in §1 bereits beschrieben, aus Bereichen mit kaltem Eis und Bereichen mit temperiertem Eis; in letzteren ist zus¨ atzlich Wasser enthalten. Unter dem Eisschild befindet sich die Lithosph¨ are, welche eine etwa 100 km dicke Schicht aus festem Gestein darstellt, die auf der viskosen Asthenosph¨ are aufschwimmt. Von der Lithosph¨ are werden im Modell allerdings nur die obersten Kilometer ber¨ ucksichtigt, da dies f¨ ur den Einfluß der thermischen Tr¨ agheit ausreichend ist. Die typische Geometrie eines polythermen Eisschildes ist in Abbildung 3.1 dargestellt, in welcher auch ein kartesisches Koordinatensystem (x, y, z) eingef¨ uhrt wird. xund y spannen hierbei die Horizontalebene auf, zist die Vertikalkoordinate antiparallel zur Richtung der Erdbeschleunigung. x, y z Atmosphere Ocean z = h(x,y,t) Asthenosphere Lithosphere Cold ice Temp. Ice CTS: z = z (x,y,t) m H(x,y,t) H (x,y,t) r z = b(x,y,t) z=b(x,y,t) r Abbildung 3.1: Skizze eines polythermen Eisschildes, in der Vertikalen stark ¨ uberh¨ oht dargestellt. Definition des in der Arbeit verwendeten kartesischen Koordinatensystems: xund ybilden die Horizontalebene, zweist nach oben. Im folgenden werden die Feldgleichungen, Randbedingungen und ¨ Ubergangsbedingungen f¨ ur dieses Problem in voller Allgemeinheit, d. h., noch ohne Vereinfachungen, entwickelt. 29 3.1 Feldgleichungen 3.1.1 Kalter Bereich Kaltes Eis bedeutet Eis mit einer Temperatur unterhalb des druckkorrigierten Schmelzpunktes. Bei Vernachl¨ assigung zus¨ atzlicher Beimischungen von Salz, Staub und Ger¨ oll kann es als viskoses w¨ armeleitendes dichtebest¨ andiges EinkomponentenFluid betrachtet werden. Die Massenbilanz lautet somit div v= 0.(3.1) Aufgrund der Annahme der Dichtebest¨ andigkeit muß der Spannungstensor Tin einen isotropen Drucktensor und einen deviatorischen Reibungsanteil aufgespalten werden, T=−p1+TR; (3.2) der Druck pist hierbei eine freie Feldgr¨ oße, wohingegen f¨ ur den Spannungsdeviator TReine Konstitutivgleichung anzugeben ist. Mit dieser Aufspaltung lautet die Impulsbilanz −grad p+ div TR+ρg=ρ˙ v.(3.3) Die noch durchzuf¨ uhrende Skalierung wird ergeben, daß der Beschleunigungsterm ρ˙ v vernachl¨ assigt werden kann, also rein Stokessches Fließen vorliegt. Es werden drei Konstitutivgleichungen ben¨ otigt: eine Spannungs-Verzerrungsgeschwindigkeits-Relation, ein Konstitutivgesetz f¨ ur die innere Energie und eines f¨ ur den W¨ armefluß q: D=EA(T0)f(σ)TRmit σ:= s1 2tr (TR)2,(3.4) ˙=c(T)˙ T, (3.5) q=−κ(T) grad T, (3.6) mit im allgemeinen temperaturabh¨ angiger spezifischer W¨ arme cund thermischer Leitf¨ ahigkeit κ. Die erste Gleichung besagt, daß die Fluidit¨ at des Eises in eine Funktion A(T0) (“Rate-Faktor”) der homologen Temperatur T0und eine Funktion f(σ) (“Kriechfunktion”) der effektiven Schubspannung σ(Wurzel der zweiten Invarianten des Spannungsdeviators IITR=1 2tr (TR)2) faktorisiert (vgl. (2.22) und Hutter [29]); die homologe Temperatur ist definiert als T0=T−TM, wobei TMdie (druckabh¨ angige) Schmelztemperatur des Eises bedeutet. Der Rate-Faktor und die Kriechfunktion sollen an dieser Stelle nicht n¨ aher spezifiziert werden; der zus¨ atzliche Faktor E (“Enhancement-Faktor”) kann auf einen Wert gr¨ oßer als Eins gesetzt werden, um 30 z. B. die leichtere Deformierbarkeit von staubhaltigem Eis gegen¨ uber normalem Eis (Paterson [62]) zu ber¨ ucksichtigen. Die zweite Gleichung verkn¨ upft ¨ Anderungen der inneren Energie ¨ uber die spezifische W¨ arme ausschließlich mit ¨ Anderungen der Temperatur, und die letzte Gleichung ist schließlich das Fouriersche W¨ armeleitgesetz. Bei Vernachl¨ assigung der Strahlungsleistung rlautet die Energiebilanz: ρ˙=−div q+ tr (TRD).(3.7) Setzt man die obigen drei Konstitutivgleichungen hier ein, so folgt hieraus eine Gleichung f¨ ur das Temperaturfeld: ρc ˙ T= div (κgrad T)+2EA(T0)f(σ)σ2.(3.8) Diese Gleichung bilanziert lokale Temperatur¨ anderungen mit Advektion (implizit in der materiellen Zeitableitung enthalten), W¨ armeleitung und dissipativer W¨ armeproduktion. 3.1.2 Temperierter Bereich Unter temperiertem Eis versteht man Eis, dessen Temperatur sich exakt auf dem druckkorrigierten Schmelzpunkt befindet, so daß dessen Temperatur nicht gesondert berechnet werden muß, sondern unmittelbar aus dem Druckfeld folgt: T=TM=T0−β∗p=T0−βp ρg,(3.9) mit T0= 0◦C sowie der Clausius-Clapeyron-Konstanten β∗(Paterson [62]) bzw. dem Clausius-Clapeyron-Gradienten β:= ρgβ∗, welcher, wie sich noch zeigen wird, dem Temperaturgradienten in temperiertem Eis entspricht. Jedoch kann temperiertes Eis einen bestimmten Anteil an Wasser enthalten; der Wassergehalt (genauer die Wasserkonzentration ω)¨ ubernimmt als thermodynamische Gr¨ oße die Rolle der Temperatur im kalten Eis. Daher muß temperiertes Eis im Gegensatz zu kaltem Eis als Mischung aus zwei verschiedenen Komponenten, n¨ amlich Eis und Wasser, angesehen werden; ρbedeutet die Dichte dieser Mischung. Es ist somit notwendig, einige Konzepte der Mischungstheorie (vgl. hierzu M¨ uller [56]) zur Anwendung zu bringen. Da allgemein angenommen wird, daß der Wassergehalt in temperierten Zonen polythermer Eisschilde mit maximal ca. 5% relativ gering ist (Hutter [31]), soll temperiertes Eis durch zwei Massenbilanzen (eine f¨ ur die Mischung als ganzes, eine f¨ ur Wasser), jedoch nur eine Impulsund Energiebilanz f¨ ur die Mischung beschrieben werden. Das Wasser wird somit als Spurenkomponente behandelt, dessen Bewegung relativ zum Baryzentrum der Mischung durch Ficksche Diffusion beschrieben wird. Alternative Konzepte, die 31 sich besser f¨ ur polytherme Gletscher mit zum Teil sehr hohem Wassergehalt eignen, wie sie in den Alpen vorkommen, verwenden statt dessen zwei Impulsbilanzen mit einer Wechselwirkungskraft vom Darcy-Typ, welche Relativbewegungen zwischen den beiden Konstituenten Eis und Wasser erm¨ oglicht (Hutter [31]). Bevor die Feldgleichungen f¨ ur temperiertes Eis formuliert werden k¨ onnen, m¨ ussen einige mischungstheoretische Gr¨ oßen eingef¨ uhrt werden. Die baryzentrische Geschwindigkeit vist definiert als v:= 1 ρ(ρivi+ρwvw).(3.10) Die Indices ibzw. wbeziehen sich auf die Komponenten Eis und Wasser, ρi/w bezeichnet die zugeh¨ orige Partialdichte. Der Wassergehalt wird als Massenkonzentration ω eingef¨ uhrt: ω:= ρw ρ.(3.11) Schließlich wird eine diffusive Wasserstromdichte jdefiniert, welche den Wasserstrom relativ zur Bewegung des Baryzentrums beschreibt: j:= ρw(vw−v) = ρω(vw−v).(3.12) Wie im Falle des kalten Eises soll auch f¨ ur temperiertes Eis Dichtebest¨ andigkeit, d. h., Konstanz von ρ, angenommen werden. Dies ist insofern problematisch, als sich die Dichten von Eis und Wasser deutlich unterscheiden (nach Paterson [62] variiert die Dichte von Gletschereis im Bereich von 830 −910 kg/m3, demgegen¨ uber steht der Wert 1000 kg/m3f¨ ur Wasser bei Normaldruck); jedoch bewegen sich aufgrund des als gering angenommenen Wassergehaltes von weniger als 5% die durch ¨ Anderungen des Wassergehaltes der Mischung Eis plus Wasser bedingten relativen Dichteschwankungen im Bereich von maximal 1% und k¨ onnen somit vernachl¨ assigt werden. Somit haben die Massenbilanz der Mischung und die Impulsbilanz der Mischung die gleiche Form wie f¨ ur kaltes Eis; sie lauten div v= 0,(3.13) −grad p+ div TR+ρg=ρ˙ v,(3.14) wobei der Spannungstensor Twiederum gem¨ aß T=−p1+TRaufgespalten wurde. Bei der Formulierung der Massenbilanz f¨ ur den Wassergehalt ist zu bedenken, daß die Partialdichte des Wassers ρwnicht konstant ist, sondern vom Wassergehalt abh¨ angt, so daß die Bilanz in der allgemeinen Form (2.5) aufgestellt werden muß. Ferner ist der Wassergehalt aufgrund der M¨ oglichkeit von Schmelzund Gefrierprozessen keine Erhaltungsgr¨ oße; es ist daher erforderlich, einen Produktionsterm Mzuzulassen: ∂ρw ∂t + div (ρwvw) = M. (3.15) 32 In ¨ aquivalenter Form kann dies als ρ˙ω=−div j+M(3.16) geschrieben werden. Wie im Falle des kalten Eises werden Konstitutivgleichungenen ben¨ otigt (vgl. auch Hutter [31]): D=EAt(ω)ft(σ)TR,(3.17) ˙=L˙ω+c(T)˙ TM,(3.18) j=−νgrad ω, (3.19) q=−κ(T) grad TM.(3.20) Die erste Gleichung, die Spannungs-Verzerrungsgeschwindigkeits-Relation, ist analog zu (3.4) f¨ ur kaltes Eis, jedoch ist die Temperaturabh¨ angigkeit des Rate-Faktors durch eine Abh¨ angigkeit vom Wassergehalt (Funktion At(ω)) ersetzt. Die zweite Gleichung besagt, daß ¨ Anderungen der inneren Energie via die latente W¨ arme Lmit ¨ Anderungen des Wassergehaltes ωund via die spezifische W¨ arme cmit ¨ Anderungen der aktuellen Schmelztemperatur zusammenh¨ angen (thermodynamisch nur approximativ, vgl. Svendsen, Greve & Hutter [69]). Die dritte Gleichung ist das bereits erw¨ ahnte Ficksche Diffusionsgesetz f¨ ur die diffusive Wasserstromdichte j, wobei νdie Diffusivit¨ at des Wassers im temperierten Eis bedeutet, und die letzte Gleichung ist schließlich das bereits bei kaltem Eis eingef¨ uhrte Fouriersche W¨ armeleitgesetz. Lund νwerden als konstant angenommen. Im folgenden soll die Energiebilanz der Mischung formuliert werden. Hierbei ist zu ber¨ ucksichtigen, daß aufgrund von (3.18) die innere Energie vom Wassergehalt ω abh¨ angt, so daß ein nichtverschwindender diffusiver Wasserstrom jeinen Fluß innerer Energie (latente W¨ arme) nach sich zieht. Somit folgt f¨ ur den gesamten W¨ armefluß qges (vgl. Gleichung (2.11)) qges =q+Lj; (3.21) der Zusatzterm Ljwurde in fr¨ uheren Arbeiten (Fowler & Larson [17], Hutter [27], Hutter [31]) nicht ber¨ ucksichtigt. Mit dieser modifizierten Form des Energieflusses folgt die Energiebilanz der Mischung zu ρ˙=−div (q+Lj) + tr (TRD).(3.22) Setzt man die Konstitutivgleichungen (3.17) – (3.20) in die Massenbilanz f¨ ur den Wassergehalt (3.16) und in die Energiebilanz (3.22) ein, so ergibt sich aus ersterer ρ˙ω=ν∇2ω+M(3.23) 33 und aus letzterer ρL ˙ω+ρc ˙ TM=Lν∇2ω+ div (κgrad TM)+2EAt(ω)ft(σ)σ2.(3.24) Konsistenz dieser beiden Gleichungen liegt vor, wenn f¨ ur die Wasserproduktion M=1 L2EAt(ω)ft(σ)σ2+ div (κgrad TM)−ρc ˙ TM(3.25) erf¨ ullt ist. Dieses Resultat ist physikalisch sinnvoll, denn es besagt, daß die W¨ armedissipation tr (TRD) im wesentlichen zum Schmelzen von Eis zu Wasser aufgebraucht, also in latente W¨ arme umgesetzt wird; anderes ist aufgrund der (n¨ aherungsweise) konstanten Temperatur auch gar nicht m¨ oglich. Abweichungen entstehen lediglich durch die kleinen Temperaturvariationen aufgrund des Clausius-Clapeyron-Gradienten. 3.1.3 Lithosph¨ are Da das Augenmerk dieser Arbeit auf der Modellierung von Eisschilden liegt, soll f¨ ur den darunter befindlichen festen Felsgrund (Lithosph¨ are) nur ein ganz einfaches Modell verwendet werden, welches lediglich die f¨ ur das Eisschild relevanten Vorg¨ ange in der Lithosph¨ are erfaßt. Hierzu wird die W¨ armeleitung und das Einsinken der Lithosph¨ are in die darunter befindliche viskose Asthenosph¨ are aufgrund des isostatischen Gleichgewichtes zwischen Eislast und Auftrieb herangezogen. Die Temperaturgleichung erh¨ alt man ganz analog zum Vorgehen bei kaltem Eis aus der Energiebilanz; sie lautet (vgl. Gl. (3.8)): ρrcr˙ T=κr∇2T. (3.26) Der Index (·)rbezieht sich jeweils auf die Lithosph¨ are, so daß ρr,crund κrderen Dichte, spezifische W¨ arme bzw. W¨ armeleitf¨ ahigkeit bedeuten. Die beiden letzteren werden im Gegensatz zu den entsprechenden Werten f¨ ur Eis als konstant angenommen. Weiterhin ist die Dissipationsleistung aufgrund m¨ oglicher Verzerrungen der Lithosph¨ are vernachl¨ assigt; nur reine W¨ armeleitungsprozesse sind ber¨ ucksichtigt. F¨ ur die Einsinktiefe ∆b(x, y, t) der Lithosph¨ are in die darunter befindliche Asthenosph¨ are aufgrund der Eislast wird zun¨ achst eine lokale Kr¨ aftebilanz zwischen Auftrieb und Eislast f¨ ur eine vertikale S¨ aule der Querschnittsfl¨ ache dA mit zugeh¨ origer Eish¨ ohe H=h−baufgestellt (aufgrund dieser Vorstellung von in der Vertikalen gegeneinander frei beweglichen S¨ aulen habe das lithosph¨ arische Geschwindigkeitsfeld keine Horizontalkomponenten): ρag∆b dA =ρgH dA, (3.27) 34 ρaist die Dichte der Asthenosph¨ are. Ist z=b0(x, y, t) die Gleichgewichtsposition der Lithosph¨ arenoberseite ohne Eislast, ergibt sich f¨ ur selbige mit Eislast, bGgw bGgw =b0−∆b=b0−ρ ρa H. (3.28) Aufgrund der Viskosit¨ at der Asthenosph¨ are stellt sich dieses Gleichgewicht nicht instantan, sondern mit einer bestimmten Verz¨ ogerungszeit τVein. F¨ ur die zeitliche Evolution der Lithosph¨ arenoberseite z=b(x, y, t) wird daher ∂b ∂t =−1 τV (b−bGgw) = −1 τV [b−(b0−ρ ρa H)] (3.29) angesetzt (Herterich [24]). Bei fester Eish¨ ohe Hentspricht dies einer exponentiellen Ann¨ aherung von ban den Gleichgewichtszustand. Unter der zus¨ atzlichen Annahme, daß jede vertikale S¨ aule der Lithosph¨ are starr ist (diese sich aber gegeneinander frei verschieben k¨ onnen) folgt, daß f¨ ur das Geschwindigkeitsfeld in der Lithosph¨ are v=∂b ∂t(x, y, t)ez(3.30) (ez: Einheitsvektor in z-Richtung) erf¨ ullt ist. 3.2 Randund ¨ Ubergangsbedingungen 3.2.1 Randbedingungen an der freien Oberfl¨ ache Wie f¨ ur jede singul¨ are Fl¨ ache l¨ aßt sich auch f¨ ur die freie Oberfl¨ ache des Eisschildes (Grenze zwischen Eis und Atmosph¨ are) eine kinematische Randbedingung formulieren. Hierzu sei angenommen, daß die freie Oberfl¨ ache durch die implizite Darstellung Fs(x, t) = 0 gegeben sei (Abbildung 3.2); die positive Seite soll mit der Atmosph¨ are, die negative mit dem Eis identifiziert werden, so daß der Normaleneinheitsvektor n= grad Fs/kgrad Fskin die Atmosph¨ are hinein zeigt. Somit muß die der Bewegung der freien Oberfl¨ ache folgende zeitliche Ableitung von Fsverschwinden: dFs dt =∂Fs ∂t +w·grad Fs= 0,(3.31) wbedeutet die Geschwindigkeit der freien Oberfl¨ ache. Die Gleichung l¨ aßt sich mit dem ¨ außeren Normaleneinheitsvektor nund dem Volumenfluß durch die freie Oberfl¨ ache hindurch a⊥ s:= (w−v−)·numschreiben zu ∂Fs ∂t +v−·grad Fs=−kgrad Fsk·a⊥ s.(3.32) 35 die negative der temperierte (untere) Bereich, so daß der Normaleneinheitsvektor n= grad Fm/kgrad Fmkin das kalte Eis hinein zeigt. Cold ice (+) Temperate ice (-) CTS: F ( ,t) = 0x m w n Abbildung 3.5: Geometrie der CTS. Zun¨ achst kann wie zuvor eine kinematische Bedingung formuliert werden. Es ergibt sich in Analogie zu (3.32) f¨ ur die freie Oberfl¨ ache ∂Fm ∂t +v·grad Fm=−kgrad Fmk·a⊥ m,(3.59) woraus mit obiger Festlegung f¨ ur Fm ∂zm ∂t +vx ∂zm ∂x +vy ∂zm ∂y −vz= 1 + ∂zm ∂x !2 + ∂zm ∂y !2  1/2 a⊥ m(3.60) folgt. In dieser Gleichung wurde der Volumenfluß durch die CTS a⊥ m:= (w−v)·neingef¨ uhrt. Diese Vorzeichenwahl bewirkt, daß a⊥ mf¨ ur Schmelzbedingungen (Str¨ omungsrichtung vom kalten in den temperierten Bereich) positiv und f¨ ur Gefrierbedingungen (Str¨ omungsrichtung vom temperierten in den kalten Bereich) negativ gez¨ ahlt wird. Aufgrund der unten abgeleiteten Kontinuit¨ at von vist es nicht n¨ otig, zwischen v+ und v−zu unterscheiden. Im Gegensatz zu der Akkumulations-Ablations-Funktion a⊥ s bei der freien Oberfl¨ ache ist a⊥ mkeine Input-Gr¨ oße, da sie im Innern des Eisschildes wirkt. a⊥ mmuß folglich vom Modell berechnet werden. Generelle Eigenschaften von Phasengrenzfl¨ achen in der Kontinuumsmechanik (vgl. Hutter [29]) sind die Kontinuit¨ at von Temperatur und Tangentialgeschwindigkeit: [[T]] = 0,[[v−(v·n)n]] = 0.(3.61) Bei der Herleitung der Massenbilanz f¨ ur temperiertes Eis wurde demonstriert, daß sich die Dichten von kaltem und temperiertem Eis um maximal 1% unterscheiden. Vernachl¨ assigt man diesen geringen Unterschied, so besagt die Massensprungbedingung 42 die Kontinuit¨ at auch der Normalgeschwindigkeit an der CTS, [[v·n]] = 0,(3.62) so daß der gesamte Geschwindigkeitsvektor kontinuierlich ist: [[v]] = 0.(3.63) Hieraus und aus der Impulssprungbedingung (2.9) ergibt sich die Kontinuit¨ at des Cauchyschen Spannungsvektors: [[T n]] = 0.(3.64) Nun soll die Massensprungbedingung f¨ ur den Wassergehalt betrachtet werden. Es ist hierbei zu ber¨ ucksichtigen, daß an der CTS Schmelzund Gefrierprozesse auftreten k¨ onnen, so daß ein Oberfl¨ achenproduktionsterm Pw mf¨ ur die Komponente Wasser eingef¨ uhrt werden muß (unten wird gezeigt, daß Pw mnur negativ oder Null sein kann, d. h., nur fl¨ achenhaftes Gefrieren, jedoch kein Schmelzen ist physikalisch m¨ oglich). Gl. (2.7) erweitert sich somit zu [[ρw(vw−w)·n]] = Pw m(3.65) oder ¨ aquivalent dazu mit der diffusiven Wasserstromdichte jgem¨ aß (3.12) (unter Ber¨ ucksichtigung der Tatsache, daß auf der positiven (kalten) Seite der CTS kein Wasser vorhanden ist, so daß die Gr¨ oßen ω+und j+gleich Null sind): −j−·n+ρa⊥ mω−=Pw m.(3.66) Dies l¨ aßt sich in eine anschaulichere Form bringen, indem statt der diffusiven Wasserstromdichte jeine totale Wasserstromdichte jtot relativ zur CTS-Bewegung w definiert wird, jtot := ρw(vw−w); hiermit schreibt sich obige Sprungbedingung als −j− tot ·n=Pw m(3.67) und besagt, daß die Normalkomponente der totalen Wasserstromdichte auf der temperierten Seite relativ zur CTS gleich der Oberfl¨ achenproduktion von Wasser ist, was unmittelbar einleuchtet. Zur Formulierung der Energiesprungbedingung gem¨ aß (2.12) muß wie f¨ ur die Ableitung von (3.22) von dem erweiterten Energiefluß in der Form (3.21) f¨ ur temperiertes Eis Gebrauch gemacht werden, so daß im kalten (positiven) Bereich qges =qund im temperierten (negativen) Bereich qges =q+Ljangesetzt wird. Unter Ber¨ ucksichtigung von (3.18), (3.63) und (3.64) erh¨ alt man q+·n−q−·n−Lj−·n=Lω−ρ(v−w)·n=−Lω−ρa⊥ m(3.68) 43 bzw. mit dem Fourierschen W¨ armeleitgesetz und der Definition von jtot κ(grad T+−grad T− M)·n+Lj− tot ·n= 0.(3.69) Da die homologe Temperatur T0=T−TMauf der kalten Seite in normaler Richtung von der CTS weg nicht ansteigen kann (sonst m¨ ußte die Temperatur die Schmelztemperatur ¨ ubersteigen), muß gelten: grad T+·n≤grad T− M·n.(3.70) Wegen (3.69) ist damit j− tot ·n≥0, so daß in der Tat f¨ ur die Wasser-Oberfl¨ achenproduktion Pw m Pw m(= −j− tot ·n)≤0 (3.71) gilt; sie kann also nicht positiv sein. Aufgrund dieser Nebenbedingung ist f¨ ur jeden Punkt der CTS zwischen drei F¨ allen zu unterscheiden, und zwar je nach dem Vorzeichen der Gr¨ oße (w−v− w)·n=a⊥ m− j−·n/(ρω−): i) (w−v− w)·n>0 (“Schmelzbedingung”): Aufgrund obiger Definition von jtot und ρw=ρω kann (3.71) nur erf¨ ullt sein, wenn ω−= 0 (3.72) gilt; in (3.71) gilt dann das Gleichheitszeichen. Einsetzen in (3.69) zeigt, daß auch grad T+·n= grad T− M·n(3.73) sein muß. Dies bedeutet, daß beim Auftreten von Schmelzbedingungen der Wassergehalt und die Normalableitung der Temperatur an der CTS stetig sind (ω+ist sowieso gleich Null, da auf der Kalteisseite der CTS nach der Definition von kaltem Eis kein Wasser vorhanden ist). ii) (w−v− w)·n<0 (“Gefrierbedingung”): In diesem Fall ist (3.71) mit ω−≥0 (3.74) vereinbar, und somit kann auch (3.70) in seiner allgemeinen Form grad T+·n≤grad T− M·n(3.75) gelten. Beim Vorliegen von Gefrierbedingungen k¨ onnen folglich der Wassergehalt und die Normalableitung der Temperatur an der CTS unstetig sein; die Unstetigkeiten 44 dieser beiden Gr¨ oßen sind ¨ uber Gl. (3.69) verkn¨ upft. iii) (w−v− w)·n= 0 (Grenzfall): Auch in diesem Fall ist (3.71) mit ω−≥0 (3.76) vertr¨ aglich, es gilt jedoch automatisch Gleichheit in (3.71). Dies in (3.69) eingesetzt ergibt grad T+·n= grad T− M·n.(3.77) Der Grenzfall ist also dadurch ausgezeichnet, daß zwar der Wassergehalt an der CTS unstetig sein kann wie bei der Gefrierbedingung, die Normalableitung der Temperatur jedoch stetig ist wie bei der Schmelzbedingung. Anschaulich kann dieses Verhalten wie folgt verstanden werden: Erreicht ein nichtverschwindender totaler Wasserstrom j− tot von der temperierten Seite her die CTS (Gefrierbedingung), so kann dieser auf der CTS ausgefrieren (negative Oberfl¨ achenproduktion von Wasser). Die hierbei freigesetzte latente W¨ arme kann dadurch abgef¨ uhrt werden, daß die Normalableitung der Temperatur auf der Kalteisseite negativer als die auf der temperierten Eisseite der CTS ist. Hiermit einher geht also ein Sprung der Normalableitung der Temperatur und (weil sich im kalten Eis kein Wasser befindet) auch des Wassergehaltes. Die umgekehrte Situation kann jedoch nicht auftreten; es ist unm¨ oglich, daß kaltes Eis auf die CTS zu fließt, zum Teil auf der CTS schmilzt (positive Oberfl¨ achenproduktion von Wasser) und so schon direkt an der CTS einen nichtverschwindenden Wasserstrom in den temperierten Bereich hinein produziert. Der Grund daf¨ ur liegt darin, daß die hierzu notwendige Schmelzw¨ arme nicht zur CTS hin transportiert werden kann, denn dazu m¨ ußte die Normalableitung der Temperatur auf der Kalteisseite positiver als die auf der temperierten Eisseite sein, was aber nicht m¨ oglich ist, da dann die Temperatur auf der Kalteisseite den Schmelzpunkt von Eis ¨ ubersteigen w¨ urde. Ein Eisfluß vom kalten in den temperierten Bereich (Schmelzbedingung) ist nur m¨ oglich ohne Oberfl¨ achenschmelzen beim Durchgang durch die CTS, so daß in diesem Fall ω−= 0 und grad T+·n= grad T− M·ngilt; in anderen Worten sind dann Wassergehalt und Temperaturgradient stetig. Es sei noch erw¨ ahnt, daß im Falle eines vernachl¨ assigbaren diffusiven Wasserstromes jim temperierten Eis, mit anderen Worten also einer sehr kleinen Wasserdiffusivit¨ at ν, die Unterscheidung zwischen Schmelzund Gefrierbedingungen einfach anhand des Vorzeichens des Eisstromes durch die CTS hindurch (a⊥ m) getroffen werden kann, da in diesem Fall v=vwgilt. a⊥ m>0 (Eisstrom vom kalten in das temperierte Eis gerichtet) entspricht dann der Schmelzbedingung, a⊥ m<0 (Eisstrom vom tempe45 rierten in das kalte Eis gerichtet) der Gefrierbedingung, und a⊥ m= 0 dem Grenzfall, was auch unmittelbar einleuchtet. Auf eine ¨ Uberpr¨ ufung des Zweiten Hauptsatzes der Thermodynamik (EntropieSprungbedingung) sei hier verzichtet. Eine solche Untersuchung findet sich in Svendsen, Greve & Hutter [69]. 46 4 Planparallele, geneigte Eisplatte 4.1 Anwendung des Modells In diesem Kapitel soll eine erste Anwendung des im vorigen Kapitel vorgestellten polythermen Eismodells gegeben werden. Wir betrachten hierzu eine zweidimensionale, in x-Richtung unendlich ausgedehnte und geneigte polytherme Eisplatte mit parallelen Seiten (“Slab”), deren Eis hangabw¨ arts fließe, wie in Abbildung 4.1 dargestellt. Numerische L¨ osungen f¨ ur eine solche Geometrie wurden bereits von Hutter, Blatter & Funk [32] und Blatter [6] konstruiert, jedoch wird hier ein etwas anderer Weg beschritten, der sogar weitgehend analytische L¨ osungen erm¨ oglicht. Atmosphere C. I. T. I. Lithosphere g z x H g z=h=H z = z (CTS) m z=b=0 Abbildung 4.1: Planparallele, geneigte, polytherme Eisplatte: Geometrie und Koordinatensystem. C. I.: kaltes Eis, T. I.: temperiertes Eis. Folgende Annahmen sollen getroffen werden: •Konstanter Neigungswinkel γund Uniformit¨ at der Prozesse in x-Richtung: (∂/∂x)(·) = 0. •Stationarit¨ at der Prozesse: (∂/∂t)(·) = 0. •Glensches Fließgesetz (vgl. Glen [19], Nye [57], Hooke [26], Paterson [62]): f(σ) = ft(σ) = σn−1(mit n= 3). •Eisfluidit¨ at unabh¨ angig von Temperatur und Wassergehalt: EA(T0) = EAt(ω)≡A= 5.3·10−24 s−1Pa−3 (Wert f¨ ur T0= 0◦C und E= 1, siehe Paterson [62]). 47 •ρ= 910 kg m−3,κ= 2.1 W m−1K−1,c= 2009 J kg−1K−1,L= 335 kJ kg−1, g= 9.81 m s−2(vgl. auch §5.5). •Vernachl¨ assigung der Druckabh¨ angigkeit des Schmelzpunktes von Eis: TM= 0◦C. •Vernachl¨ assigung der Wasserdiffusion: ν= 0 ⇒j=0. •Vernachl¨ assigung von lithosph¨ arischen Einfl¨ ussen: kein Einsinken in die Asthenosph¨ are, keine Berechnung von Temperatur und basaler Schmelzrate. Mit diesen Annahmen lauten die Modellgleichungen wie folgt: Massenbilanz, kalter und temperierter Bereich (aus Gln. (3.1), (3.13)): dvz dz = 0.(4.1) Impulsbilanz, kalter und temperierter Bereich (aus Gln. (3.3), (3.14) mit vernachl¨ assigter Beschleunigung): dσxz dz +ρg sin γ= 0,(4.2) −dp dz +dσR z dz −ρg cos γ= 0.(4.3) Energiebilanz, kalter Bereich (aus Gl. (3.8)): ρcvz dT dz =κd2T dz2+ 2A σ4.(4.4) Energiebilanz, temperierter Bereich, bzw. Massenbilanz f¨ ur den Wassergehalt (aus Gln. (3.23), (3.24)): ρvz dω dz = 2A Lσ4.(4.5) Spannungs-Verzerrungsgeschwindigkeits-Relation, kalter und temperierter Bereich (aus Gln. (3.4), (3.17)): σR x= 0,(4.6) σR z= 0,(4.7) dvx dz = 2A σ2σxz; (4.8) hieraus ergibt sich mit der Definition der effektiven Schubspannung σ:= qtr (TR)2/2 diese zu σ=σxz. 48 Randbedingungen, kalte freie Oberfl¨ ache (aus Gln. (3.33), (3.34), (3.35)): vz=−a⊥ s,(4.9) σ=σxz = 0,(4.10) −p+σR z=−p= 0,(4.11) T=Ts.(4.12) Randbedingungen, temperierte Basis: Aufgrund der Gleichungen (4.1) und (4.9) ist die Vertikalgeschwindigkeit vzauf der gesamten H¨ ohe des Slabs gleich der negativen Akkumulations-Ablations-Funktion a⊥ s. Insbesondere nimmt vzalso auch an der Basis diesen Wert an, die kinematische Bedingung (3.50) wird folglich gar nicht ben¨ otigt. Man k¨ onnte aus dieser Gleichung den Wasser-Massenstrom ˙mw bin den Boden hinein berechnen, der erforderlich w¨ are, um dieses basale vzaufrechtzuerhalten; dies ist jedoch wenig interessant und wird daher nicht durchgef¨ uhrt. Randbedingung (3.51), welche die Normalkomponente der diffusiven Wasserstromdichte bestimmt, ist f¨ ur das Slab-Problem auch nicht erforderlich, da die Wasserdiffusion sowieso vernachl¨ assigt wird. Somit verbleibt das Gleitgesetz (3.52), um eine Randbedingung f¨ ur die basale Tangentialgeschwindigkeit vx,b zu erhalten. Der Einfachheit halber soll jedoch auf die explizite Angabe eines Gleitgesetzes verzichtet werden; statt dessen wird vx,b selbst vorgeschrieben. Die einzige Auswirkung des Wertes von vx,b auf die Resultate besteht ohnehin darin, daß sich vxals Funktion der H¨ ohe z um eine additive Konstante ¨ andert; das Verhalten der Temperatur und des Wassergehaltes wird nicht beeinflußt. ¨ Ubergangsbedingungen, CTS (aus Gl. (3.61), (3.63), (3.64), (3.68), (3.70)): T+=T−,(4.13) v+ x=v− x, v+ z=v− z,(4.14) p+=p−, σ+ (xz)=σ− (xz),(4.15) κdT+ dz =Lω−ρa⊥ mmit dT+ dz ≤0.(4.16) Die Nebenbedingung in letzterer Gleichung bewirkt, daß zwei verschiedene F¨ alle unterschieden werden m¨ ussen (vgl. Diskussion in §3.2.5): •a⊥ m>0 (“Schmelzbedingung”, Eisfluß vom kalten in den temperierten Bereich): dT+/dz = 0, ω−= 0. •a⊥ m<0 (“Gefrierbedingung”, Eisfluß vom temperierten in den kalten Bereich): Gleichung (4.16) in ihrer nichttrivialen Form, d. h., dT+/dz kann strikt negativ 49 und ω−strikt positiv sein; in diesem Fall ist eine zus¨ atzliche Randbedingung f¨ ur den basalen Wassergehalt erforderlich. Der Grenzfall a⊥ m= 0 wird nicht betrachtet, da er keinen station¨ aren Zustand zul¨ aßt. Wegen (4.21) w¨ are n¨ amlich vz≡0, somit m¨ ußte die linke Seite in (4.5) verschwinden, gleichzeitig aber deren rechte Seite wegen (4.19) gr¨ oßer als Null sein, was einen Widerspruch darstellt. 4.2 Integration der Slab-Gleichungen Die oben hergeleiteten Gleichungen, welche aus der Anwendung des polythermen Eismodells aus Kapitel 3 auf das spezielle Slab-Problem hervorgegangen sind, k¨ onnen fast vollst¨ andig analytisch integriert werden. Lediglich f¨ ur die Position z=zmder CTS verbleibt eine nichtlineare algebraische Gleichung, welche aber leicht mit Hilfe des Newton-Verfahrens zur Nullstellenbestimmung gel¨ ost werden kann. Dies soll im folgenden durchgef¨ uhrt werden. Integration von (4.2) und (4.3) ergibt unter Ber¨ ucksichtigung von (4.7), (4.10), (4.11) und (4.15) p(z) = ρg cos γ(H−z),(4.17) σxz(z) = ρg sin γ(H−z) (4.18) und somit σ=σxz =ρg sin γ(H−z).(4.19) Der Druck verh¨ alt sich offensichtlich rein hydrostatisch. Das Geschwindigkeitsfeld ergibt sich mit Hilfe dieses Ergebnisses aus (4.1), (4.8), (4.9), (4.14) sowie vorgeschriebenem vx,b zu vx(z) = A 2(ρg sin γ)3[H4−(H−z)4] + vx,b,(4.20) vz(z) = const = −a⊥ s=−a⊥ m.(4.21) Die horizontale Geschwindigkeit vxsteigt monoton von ihrem minimalen Wert vx,b an der Basis zu einem maximalen Wert an der freien Oberfl¨ ache an, wie man es von einer Scherstr¨ omung mit freier Oberfl¨ ache erwartet. Die vertikale Geschwindigkeit vz ist, wie oben bereits bemerkt wurde, ¨ uber die H¨ ohe hinweg konstant. Sie bilanziert an der freien Oberfl¨ ache die Akkumulations-Ablations-Funktion a⊥ s, und an der CTS entspricht sie dem (negativen) Volumenfluß a⊥ mdurch die CTS hindurch. Die L¨ osung der Gleichungen (4.4) und (4.5) f¨ ur den Temperaturund Wassergehalt im kalten bzw. temperierten Bereich sowie die damit verbundene Bestimmung 50 der CTS-Position erfordert einen ziemlich großen Rechenaufwand. Zun¨ achst wird Gleichung (4.19) f¨ ur σund (4.21) f¨ ur vzeingesetzt; dies ergibt κd2T dz2+ρca⊥ s dT dz =−2A(ρg sin γ)4(H−z)4(4.22) und ρa⊥ s dω dz =−2A L(ρg sin γ)4(H−z)4.(4.23) Zwecks einfacherer Rechnung wird nun die vertikale Koordinate zmit der Transformation z=Hζ auf das Intervall [0,1] abgebildet. Man erh¨ alt so Dd2T dζ2+MdT dζ =−K(1 −ζ)4(4.24) und Mdω dζ =−Kt(1 −ζ)4(4.25) mit den Abk¨ urzungen D=κ ρc, M=Ha⊥ s, K=2A ρc H6(ρg sin γ)4, Kt=2A ρLH6(ρg sin γ)4.(4.26) Zuerst soll die Temperaturgleichung (4.24) betrachtet werden. Zun¨ achst muß die homogene Gleichung gel¨ ost werden; mit dem Ansatz T∼exp(λζ) ergibt sich Dλ2+Mλ = 0 ⇒λ1= 0, λ2=−M D.(4.27) Ein Partikularintegral der inhomogenen Gleichung muß in der Form Tp=a1ζ+a2ζ2+a3ζ3+a4ζ4+a5ζ5(4.28) existieren; die Koeffizienten a1bis a5erh¨ alt man durch Einsetzen hiervon in die Temperaturgleichung (4.24) und nachfolgenden Koeffizientenvergleich der Terme von der Ordnung [ζ0] bis [ζ4]. Dies liefert ein lineares Gleichungssystem: ζ0: 2Da2+Ma1=−K, ζ1: 6Da3+ 2Ma2= 4K, ζ2: 12Da4+ 3Ma3=−6K, ζ3: 20Da5+ 4Ma4= 4K, ζ4: 5Ma5=−K; (4.29) 51 Abbildung 4.3: Horizontalgeschwindigkeit vx(Plot a) und Temperatur Tbzw. Wassergehalt ω (Plot b) f¨ ur den Slab mit Gefrierbedingung an der CTS. H= 200 m, γ= 4◦,Ts=−10◦C, a⊥ s= −0.2 m a−1,vx,b = 5 m a−1. 58 5 Skalierung und Flacheisannahme Das in Kapitel 3 vorgestellte polytherme Eismodell ist f¨ ur einen numerischen oder gar analytischen Zugang im allgemeinen noch zu kompliziert. Lediglich sehr einfache Probleme mit hochgradiger Symmetrie (z. B. geneigte planparallele Eisplatte, siehe voriges Kapitel) k¨ onnen damit behandelt werden. Es ist also notwendig, die Gleichungen weiter zu vereinfachen. Hierzu werden sie zun¨ achst einer Skalierung unterworfen, wobei die physikalischen Gr¨ oßen durch Wahl von zugeh¨ origen typischen Werten entdimensioniert werden; dies erm¨ oglicht in einem weiteren Schritt systematische Vernachl¨ assigungen von Termen aufgrund der Gr¨ oßenordnung der entstehenden dimensionslosen Produkte. In unserem konkreten Fall beruhen die Vernachl¨ assigungen im wesentlichen auf der Annahme, daß die im Eisschild auftretenden typischen H¨ ohen wesentlich kleiner als die typischen L¨ angendimensionen sind, in anderen Worten auf der Annahme eines hinreichend flachen Eisschildes (“Flacheisannahme” bzw. “Shallow-ice-approximation”, im folgenden als “SIA” bezeichnet; vgl. Hutter [29], Morland [55], Blatter [6]). 5.1 Einf¨ uhrung der Skalierung Die verschiedenen physikalischen Gr¨ oßen sollen wie folgt skaliert werden: (x, y) = [L] (˜x, ˜y), z= [H] ˜z, (vx, vy) = [VL] (˜vx,˜vy), vz= [VH] ˜vzmit ε:= [H] [L]=[VH] [VL], t=[L] [VL]˜ t, (T, T0) = [∆T] (˜ θ, ˜ θ0) (T, T0in ◦C), ω= [ω] ˜ω, p=ρg[H] ˜p, (σxz, σyz, σ) = ερg[H] (˜σxz,˜σyz,˜σ), (σR x, σR y, σR z, σxy) = ε2ρg[H] (˜σR x,˜σR y,˜σR z,˜σxy), (h, zm, b, br) = [H] (˜ h, ˜zm,˜ b,˜ br), (a⊥ s, a⊥ m) = [VH] (˜a⊥ s,˜a⊥ m), (A(T0), At(ω)) = [A] ( ˜ A(˜ θ0),˜ At(˜ω)), (f(σ), ft(σ)) = [f] ( ˜ f(˜σ),˜ ft(˜σ)), 59 Q⊥ geoth = [Q⊥ geoth]˜ Q⊥ geoth, κ(T) = [κ] ˜κ(˜ θ), κr= [κr] ˜κr, c(T) = [c] ˜c(˜ θ), cr= [cr] ˜cr, C(t⊥, . . .) = [C]˜ C(˜ t⊥, . . .), Ct(t⊥, . . .) = [Ct]˜ Ct(˜ t⊥, . . .), Pw b=ρ[ω][VH]˜ Pw b, ˙mw b=ρ[ω][VH]˜ ˙mw b, Pw m=ρ[ω][VH]˜ Pw m.(5.1) Gr¨ oßen in eckigen Klammern bedeuten typische Werte f¨ ur die jeweiligen Variablen, und mit einer Tilde versehene Variablen sind dimensionslos. Wenn nicht ausdr¨ ucklich anders erw¨ ahnt, werden im folgenden nur noch dimensionslose Gr¨ oßen verwendet; daher werden die Tilden der Einfachheit halber wieder weggelassen. Die Skalierung soll so gew¨ ahlt sein, daß die dimensionslosen Gr¨ oßen allesamt von der Ordnung O(1) sind. Dies erscheint bei der zun¨ achst etwas fragw¨ urdigen Skalierung der verschiedenen Spannungen zweifelhaft, jedoch wird sich weiter unten zeigen, daß (ausgehend von der sicherlich vern¨ unftigen Skalierung des Druckes pmit dem ¨ Uberlagerungsdruck ρg[H] und der Annahme, daß von den sechs deviatorischen Spannungen die Schubspannungen σxz und σyz die dominierenden sind) genau diese Skalierung die richtige ist, d. h., Konsistenz in den Ordnungen von linken und rechten Seiten der Terme der skalierten Impulsbilanzen und der skalierten SpannungsVerzerrungsgeschwindigkeits-Relationen herstellt. Diese Spezifizierung ist ad¨ aquat f¨ ur die hier behandelten Eisschilde, welche auf festem Grund aufliegen, so daß die Schubspannungen σxz,σyz am Boden ¨ ubertragen werden k¨ onnen. Bei Schelfeisen, d. h., aufschwimmenden Eismassen, sind dagegen gerade diese Schubspannungen nicht relevant, da Wasser als n¨ aherungsweise ideales Fluid keine Schubspannungen aufrecht erhalten kann. Somit ist dies der Punkt in der Herleitung des Modells, welcher die Aufspaltung zwischen Eisschilden und Schelfeisen darstellt. Die oben eingef¨ uhrten 13 typischen Gr¨ oßen [L], [H], [VL], [∆T], [ω], [A][f] (tritt nur als Produkt auf), [Q⊥ geoth], [κ], [κr], [c], [cr], [C], [Ct] ([VH] ist nicht unabh¨ angig und wird daher nicht aufgef¨ uhrt) bilden zusammen mit den acht physikalischen Konstanten ρ,g,L,β,ν,ρr,ρa,τVeinen Satz von 21 Gr¨ oßen, deren Dimensionen aus den vier Grundeinheiten Meter, Kilogramm, Sekunde und Kelvin bestehen. Die zugeh¨ orige 60 Dimensionsmatrix hat offensichtlich den Rang vier: [L]ρ[VL] [∆T]. . . m 1 −3 1 0 . . . kg 0 1 0 0 . . . s 0 0 −1 0 . . . K 0 0 0 1 . . . Nach den Regeln der Dimensionsanalyse (Hutter [30]) existieren also 21 −4 = 17 unabh¨ angige dimensionslose Produkte. Ein vollst¨ andiger Satz solcher Produkte ist ε=[H] [L]=[VH] [VL], F=[VL]2 g[L], D=[κ] ρ[c][H][VH], α=g[H] [c][∆T], K=ρg[H]3[A][f] [L][VL], αt=g[H] L[ω], Dt=ν ρ[H][VH], F=ρg[H]2[C] [L][VL], Ft=ρg[H]2[Ct] [L][VL], B=β[H] [∆T], Dr=[κr] ρr[cr][H][VH], Tr=τV[VH] [H], Nr=[H][Qgeoth] [κr][∆T], [ω],[κr] [κ],ρr ρ,ρa ρ; (5.2) εist das Aspektverh¨ altnis, Fdie Froude-Zahl, Ddie W¨ armediffusionszahlzahl, αdas Verh¨ altnis potentielle Energie zu innere Energie in kaltem Eis, Kdie Fluidit¨ atszahl, 61 αtdas Verh¨ altnis potentielle Energie zu innere Energie in temperiertem Eis, Dtdie Wasserdiffusionszahl, Fbzw. Ftdie Gleitzahl f¨ ur kalte bzw. temperierte Basis, B die Clausius-Clapeyron-Zahl, Drdie W¨ armediffusionszahl der Lithosph¨ are, Trdie Zeitverz¨ ogerungszahl f¨ ur das Einsinken der Lithosph¨ are in die Asthenosph¨ are, und schließlich Nrdie geotherme W¨ armezahl in der Lithosph¨ are. Die vier letzten Kombinationen erhalten keine speziellen Bezeichnungen. 5.2 Skalierung und SIA f¨ ur die Modellgleichungen Im folgenden geht es darum, alle Feldgleichungen, Randbedingungen und ¨ Ubergangsbedingungen des Kapitels 3 obiger Skalierung zu unterziehen und sie anschließend mit der SIA (Flacheisannahme) zu vereinfachen. Wie bereits erw¨ ahnt besteht diese SIA in der Annahme eines hinreichend flachen Eisschildes, dessen typische L¨ angenabmessungen wesentlich gr¨ oßer als dessen typische H¨ ohenabmessungen sind, was in der Natur meist gut erf¨ ullt ist (typische Werte liegen im Bereich von ε= 10−3). Formelhaft ausgedr¨ uckt bedeutet dies, daß der Limes ε→0 (5.3) betrachtet wird, also alle Terme, die von der Ordnung O(εp) mit p≥1 sind, in den Gleichungen zu vernachl¨ assigen sind. Weiterhin l¨ aßt sich leicht absch¨ atzen, daß die Froude-Zahl Fbei realen Eisschilden in der Gr¨ oßenordnung F= 10−15 und weniger liegt, es also sinnvoll ist, gleichzeitig den Limes F ε→0 (5.4) zu betrachten. 5.2.1 Kalter Bereich Die skalierte Massenbilanz ergibt sich durch Einsetzen der Skalierung in die urspr¨ ungliche Massenbilanz (3.1) zu ∂vx ∂x +∂vy ∂y +∂vz ∂z = 0 (5.5) (wie bereits erw¨ ahnt werden die Tilden, welche in Abschnitt 5.1 dimensionslose Gr¨ oßen anzeigen, im folgenden weggelassen, da nur noch mit solchen gerechnet wird), und kann durch die SIA nicht weiter vereinfacht werden. Aus der Impulsbilanz (3.3) folgt F ε dvx dt =−∂p ∂x +ε2∂σR x ∂x +ε2∂σxy ∂y +∂σxz ∂z ,(5.6) 62 F ε dvy dt =ε2∂σxy ∂x −∂p ∂y +ε2∂σR y ∂y +∂σyz ∂z ,(5.7) Fεdvz dt =ε2∂σxz ∂x +ε2∂σyz ∂y −∂p ∂z +ε2∂σR z ∂z −1.(5.8) Anwendung der SIA ergibt weiter −∂p ∂x +∂σxz ∂z = 0,(5.9) −∂p ∂y +∂σyz ∂z = 0,(5.10) −∂p ∂z = 1; (5.11) dies gibt uns ein erstes, sehr eindrucksvolles Beispiel, wie sehr das Modell des Kapitels 3 durch die SIA vereinfacht wird. Das Ergebnis zeigt ferner, daß die Skalierung der Schubspannungen σxz und σyz mit ερg[H] (vgl. Abschnitt 5.1) richtig war, denn ohne das zus¨ atzliche εim Vergleich zur Druckskalierung w¨ aren die Gleichungen (5.9) und (5.10) nicht konsistent in den Ordnungen der Summanden, sondern ein Term der Ordnung O(ε) und ein Term der Ordnung O(1) m¨ ußten sich zu Null aufsummieren, was nicht m¨ oglich ist. Die Energiebilanz (3.8) ergibt mit ausgeschriebenem Advektionsterm ∂θ ∂t +vx ∂θ ∂x +vy ∂θ ∂y +vz ∂θ ∂z =D c(ε2∂ ∂x κ∂θ ∂x! +ε2∂ ∂y κ∂θ ∂y!+∂ ∂z κ∂θ ∂z !)+ 2α cKEA(θ0)f(σ)σ2; (5.12) mit der SIA wird hieraus ∂θ ∂t +vx ∂θ ∂x +vy ∂θ ∂y +vz ∂θ ∂z =D c ∂ ∂z κ∂θ ∂z !+ 2α cKEA(θ0)f(σ)σ2.(5.13) Die horizontale W¨ armeleitung kann also im Rahmen des SIA-Limes vernachl¨ assigt werden, jedoch muß sowohl die horizontale als auch die vertikale Advektion ber¨ ucksichtigt werden. Von den Materialgesetzen (3.4) – (3.6), die bereits in die Energiebilanz eingearbeitet sind, wird im folgenden nur noch die Spannungs-VerzerrungsgeschwindigkeitsRelation (3.4) ben¨ otigt: ∂vx ∂x =KEA(θ0)f(σ)σR x,(5.14) ∂vy ∂y =KEA(θ0)f(σ)σR y,(5.15) ∂vz ∂z =KEA(θ0)f(σ)σR z,(5.16) 63 ∂vx ∂y +∂vy ∂x = 2KEA(θ0)f(σ)σxy,(5.17) ∂vx ∂z +ε2∂vz ∂x = 2KEA(θ0)f(σ)σxz,(5.18) ∂vy ∂z +ε2∂vz ∂y = 2KEA(θ0)f(σ)σyz.(5.19) Mit der SIA vereinfacht sich dies zu ∂vx ∂x =KEA(θ0)f(σ)σR x,(5.20) ∂vy ∂y =KEA(θ0)f(σ)σR y,(5.21) ∂vz ∂z =KEA(θ0)f(σ)σR z,(5.22) ∂vx ∂y +∂vy ∂x = 2KEA(θ0)f(σ)σxy,(5.23) ∂vx ∂z = 2KEA(θ0)f(σ)σxz,(5.24) ∂vy ∂z = 2KEA(θ0)f(σ)σyz.(5.25) Hieraus wird ersichtlich, daß auch die Skalierung der normalen Reibungsspannungen σR x,σR yund σR zsowie der Schubspannung σxy gem¨ aß Abschnitt 5.1 mit ε2ρg[H] vern¨ unftig war, also ordnungskonsistente Terme ergibt. Diese Spannungen sind folglich f¨ ur flache Eisschilde sehr klein, was auch intuitiv unmittelbar einleuchtet. Da im folgenden nur die letzten beiden dieser sechs Gleichungen ben¨ otigt werden (n¨ amlich zur Berechnung der Horizontalgeschwindigkeiten), brauchen sie nicht weiter betrachtet zu werden. Schließlich erh¨ alt man f¨ ur die in der Energiebilanz und im Materialgesetz ben¨ otigte effektive Schubspannung σ=v u u tσ2 xz +σ2 yz +ε2 (σR x)2 2+(σR y)2 2+(σR z)2 2+σ2 xy!; (5.26) im SIA-Limes ergibt dies σ=qσ2 xz +σ2 yz.(5.27) 5.2.2 Temperierter Bereich Wie bereits festgestellt, haben die Massenbilanz der Mischung Eis plus Wasser und die Impulsbilanz der Mischung die gleiche Form wie die entsprechenden Bilanzen f¨ ur kaltes Eis (vgl. (3.1), (3.3), (3.13), (3.14)); daher ¨ ubertragen sich auch die skalierten und die SIA-Versionen hiervon aus dem vorigen Abschnitt (Gln. (5.5) – (5.11)). Sie werden hier nicht noch einmal explizit aufgef¨ uhrt. 64 F¨ ur die Temperatur gilt gem¨ aß (3.9) θ=θM=−Bp. (5.28) Aus der Energiebilanz der Mischung (3.24) bzw. der Massenbilanz f¨ ur den Wassergehalt (3.23) (diese beiden Gleichungen sind identisch) erh¨ alt man ∂ω ∂t +vx ∂ω ∂x +vy ∂ω ∂y +vz ∂ω ∂z +cαt α ∂θM ∂t +vx ∂θM ∂x +vy ∂θM ∂y +vz ∂θM ∂z ! =Dt ε2∂2ω ∂x2+ε2∂2ω ∂y2+∂2ω ∂z2! +Dαt α(ε2∂ ∂x κ∂θM ∂x !+ε2∂ ∂y κ∂θM ∂y !+∂ ∂z κ∂θM ∂z !) +2αtKEAt(ω)ft(σ)σ2,(5.29) was sich durch die SIA vereinfacht zu ∂ω ∂t +vx ∂ω ∂x +vy ∂ω ∂y +vz ∂ω ∂z +cαt α ∂θM ∂t +vx ∂θM ∂x +vy ∂θM ∂y +vz ∂θM ∂z ! =Dt ∂2ω ∂z2+Dαt α ∂ ∂z κ∂θM ∂z !+ 2αtKEAt(ω)ft(σ)σ2.(5.30) Wie bei der Energiebilanz f¨ ur kaltes Eis f¨ allt die horizontale Diffusion durch die SIA heraus. Die Spannungs-Verzerrungsgeschwindigkeits-Relation (3.17) hat im wesentlichen die gleiche Form wie diejenige f¨ ur kaltes Eis (3.4). Daher gelten die Gleichungen (5.14) – (5.25) auch f¨ ur temperiertes Eis; es ist lediglich A(θ0) durch At(ω) und f(σ) durch ft(σ) zu ersetzen. Ebenso bleibt das Ergebnis f¨ ur die effektive Schubspannung unver¨ andert. 5.2.3 Lithosph¨ are Die Skalierung der Energiebilanz (3.26) ergibt ∂θ ∂t +vx ∂θ ∂x +vy ∂θ ∂y +vz ∂θ ∂z =Drκr cr ε2∂2θ ∂x2+ε2∂2θ ∂y2+∂2θ ∂z2!,(5.31) bzw. im SIA-Limes ∂θ ∂t +vx ∂θ ∂x +vy ∂θ ∂y +vz ∂θ ∂z =Drκr cr ∂2θ ∂z2.(5.32) Die zeitliche Evolution der Lithosph¨ arenoberseite folgt aus (3.29) zu ∂b ∂t =−1 Tr [b−(b0−ρ ρa H)]; (5.33) 65 dies kann nicht weiter vereinfacht werden. Entsprechend gilt f¨ ur die Geschwindigkeit der Lithosph¨ are gem¨ aß (3.30) vx= 0,(5.34) vy= 0,(5.35) vz=∂b ∂t(x, y, t) (5.36) mit ∂b/∂t aus obiger Gleichung. 5.2.4 Randbedingungen an der freien Oberfl¨ ache Zun¨ achst soll die kinematische Bedingung (3.33) skaliert werden. Dies ergibt3 ∂h ∂t +vx ∂h ∂x +vy ∂h ∂y −vz= 1 + ε2 ∂h ∂x!2 +ε2 ∂h ∂y !2  1/2 a⊥ s; (5.37) mit der SIA verbleibt ∂h ∂t +vx ∂h ∂x +vy ∂h ∂y −vz=a⊥ s.(5.38) F¨ ur das weitere Vorgehen werden die Terme grad Fsund kgrad Fskben¨ otigt. Man erh¨ alt hierf¨ ur grad Fs= −ε∂h ∂x,−ε∂h ∂y ,1!t (5.39) und kgrad Fsk= 1 + ε2 ∂h ∂x!2 +ε2 ∂h ∂y !2  1/2 ; (5.40) im SIA-Limes ergibt dies grad Fs= (0,0,1)t,kgrad Fsk= 1.(5.41) Mit dem ¨ außeren Normalenvektor n= grad Fs/kgrad Fskfolgt aus der Impulssprungbedingung (3.34) −(−p+ε2σR x)∂h ∂x −ε2σxy ∂h ∂y +σxz = 0,(5.42) −ε2σxy ∂h ∂x −(−p+ε2σR y)∂h ∂y +σyz = 0,(5.43) −ε2σxz ∂h ∂x −ε2σyz ∂h ∂y −p+ε2σR z= 0.(5.44) 3Die hochgestellten Minuszeichen (·)−zur Kennzeichnung der Eis-Seite bei den Randbedingungen an der freien Oberfl¨ ache sowie den ¨ Ubergangsbedingungen an der Eisbasis (vgl. Kap. 3) werden im folgenden weggelassen. Nicht gekennzeichnete Gr¨ oßen beziehen sich daher stets auf die Eis-Seite. 66 Anwendung der SIA vereinfacht zun¨ achst die letzte dieser Gleichungen zu p= 0,(5.45) und hiermit ergeben die beiden ersten Gleichungen σxz = 0, σyz = 0.(5.46) Es verbleibt im Falle einer kalten freien Oberfl¨ ache die Randbedingung f¨ ur die Oberfl¨ achentemperatur (3.35), welche skaliert einfach θ(x, y, t) = θs(x, y, t) (5.47) ergibt; dies kann offensichtlich nicht weiter vereinfacht werden. Im Falle einer temperierten freien Oberfl¨ ache ist in analoger Weise der Wassergehalt ωoder dessen Normalableitung vorzuschreiben. 5.2.5 ¨ Ubergangsbedingungen an der kalten Eisbasis Die kinematische Bedingung (3.38) wird durch die Skalierung und SIA nicht ver¨ andert, sie lautet ∂b ∂t +vx ∂b ∂x +vy ∂b ∂y −vz= 0.(5.48) Analog zur freien Oberfl¨ ache seien nun die Terme grad Fbund kgrad Fbkberechnet: grad Fb= ε∂b ∂x, ε∂b ∂y,−1!t (5.49) und kgrad Fbk= 1 + ε2 ∂b ∂x!2 +ε2 ∂b ∂y!2  1/2 ; (5.50) im SIA-Limes erh¨ alt man hieraus grad Fb= (0,0,−1)t,kgrad Fbk= 1.(5.51) Der ¨ außere Normalenvektor folgt wie zuvor zu n= grad Fb/kgrad Fbk. Mit T n =ρg[H] kgrad Fbk            εσx ∂b ∂x +ε3σxy ∂b ∂y −εσxz ε3σxy ∂b ∂x +εσy ∂b ∂y −εσyz ε2σxz ∂b ∂x +ε2σyz ∂b ∂y −σz            (5.52) 67 5.3.1 Berechnung der Spannungen Aus Gleichung (5.11) zusammen mit der Randbedingung (5.45) und der ¨ Ubergangsbedingung (5.88) erh¨ alt man f¨ ur den Druck p(x, y, z, t) = h(x, y, t)−z; (5.95) er unterliegt also einer rein hydrostatischen Verteilung. Unter Verwendung dieses Ergebnisses folgt aus (5.9), (5.10), (5.46) und (5.88) σxz =−∂h ∂x(h−z),(5.96) σyz =−∂h ∂y (h−z).(5.97) Diese Schubspannungen ergeben sich also als Produkt aus Oberfl¨ achenneigung und ¨ Uberlagerungsdruck. F¨ ur die effektive Schubspannung ergibt sich hiermit gem¨ aß (5.27) σ= (h−z)v u u t ∂h ∂x!2 + ∂h ∂y !2 .(5.98) Man merke sich, daß die obigen Ergebnisse gleichermaßen im kalten und temperierten Bereich des Eisschildes gelten. Die verbleibenden, sehr kleinen Spannungen σR x,σR y,σR zund σxy k¨ onnen nicht direkt berechnet werden, sie werden jedoch auch im weiteren nicht ben¨ otigt. Man kann sie im Anschluß an die Berechnung des Geschwindigkeitsfeldes aus den ersten vier Gleichungen (5.20) – (5.23) der Spannungs-Verzerrungsgeschwindigkeits-Relation ermitteln. 5.3.2 Berechnung der Geschwindigkeit Die Horizontalgeschwindigkeiten vxund vylassen sich durch Einsetzen der obigen Ergebnisse f¨ ur σxz und σyz in die beiden letzten Gleichungen (5.24), (5.25) der Spannungs-Verzerrungsgeschwindigkeits-Relation und anschließende Integration erhalten: vx=vx,b −2K∂h ∂x Zz bEA(t)(·)f(t)(σ)(h−z0)dz0,(5.99) vy=vy,b −2K∂h ∂y Zz bEA(t)(·)f(t)(σ)(h−z0)dz0.(5.100) Hierbei ist A(t)(·) =    A(θ0) f¨ ur z > zm(kalter Bereich), At(ω) f¨ ur z < zm(temperierter Bereich) (5.101) 74 und f(t)(σ) =    f(σ) f¨ ur z > zm(kalter Bereich), ft(σ) f¨ ur z < zm(temperierter Bereich).(5.102) Die Integrationskonstanten vx,b und vy,b, welche die entsprechenden basalen Geschwindigkeiten darstellen, folgen aus dem Gleitgesetz (5.57), (5.58) unter Verwendung der Ergebnisse f¨ ur σxz und σyz zu vx,b = (vsl)x=−F(t)C(t) ∂h ∂x(h−b),(5.103) vy,b = (vsl)y=−F(t)C(t) ∂h ∂y (h−b).(5.104) Hierbei wurde ausgenutzt, daß die Lithosph¨ arengeschwindigkeit keine xund y-Komponenten hat (vgl. (5.34), (5.35)), so daß gem¨ aß (5.63), (5.64) die entsprechenden Gleitgeschwindigkeitskomponenten mit den Eisgeschwindigkeitskomponenten identisch sind. Weiterhin wurde F(t)C(t)=  FCf¨ ur kalte Basis, FtCtf¨ ur temperierte Basis (5.105) gesetzt. Mit den obigen Ergebnissen f¨ ur die Horizontalgeschwindigkeiten kann man aus der Massenbilanz (5.5) auch die Vertikalgeschwindigkeit berechnen: vz=−Zz b ∂vx ∂x +∂vy ∂y !dz0+vz,b.(5.106) Die basale Vertikalgeschwindigkeit vz,b kann gem¨ aß der kinematischen Bedingung (5.48) bzw. (5.70) durch vx,b und vy,b ausgedr¨ uckt werden. Zusammenfassend gilt also f¨ ur die Geschwindigkeiten: vx=−F(t)C(t) ∂h ∂x(h−b)−2K∂h ∂x Zz bEA(t)(·)f(t)(σ)(h−z0)dz0,(5.107) vy=−F(t)C(t) ∂h ∂y (h−b)−2K∂h ∂y Zz bEA(t)(·)f(t)(σ)(h−z0)dz0,(5.108) vz=−Zz b ∂vx ∂x +∂vy ∂y !dz0+∂b ∂t +vx,b ∂b ∂x +vy,b ∂b ∂y −[ω] ˙mw b.(5.109) Bei bekanntem Temperaturund Wassergehaltsfeld k¨ onnen alle Integrale leicht numerisch berechnet werden. Nimmt man die Eisfluidit¨ at als temperaturbzw. wassergehaltsunabh¨ angig an (A(t)= const), so treten Abh¨ angigkeiten von Temperatur und Wassergehalt gar nicht mehr auf, und die Integrationen k¨ onnen sogar analytisch durchgef¨ uhrt werden. 75 Abschließend sei noch bemerkt, daß offensichtlich (vx, vy)∝ −(∂h/∂x, ∂h/∂y) = −grad(x,y)hgilt, d. h., der Vektor der Horizontalgeschwindigkeit zeigt in jedem Punkt des Eisschildes in Richtung des st¨ arksten Oberfl¨ achengef¨ alles, insbesondere also entlang eines Vertikalprofils immer in die gleiche Richtung (vgl. Hutter [29]). Dieses ¨ uberraschende Ergebnis erm¨ oglicht es, durch Bestimmung der Richtung der Eisflußgeschwindigkeit in Eisbohrkernen die G¨ ultigkeit der SIA in einem gegebenen realen Eisschild zu ¨ uberpr¨ ufen. 5.3.3 Evolution der freien Oberfl¨ ache Eine Gleichung, die die zeitliche Evolution der freien Oberfl¨ ache, mit anderen Worten also die Dicken¨ anderung des Eisschildes mit der Zeit beschreibt, kann aus der Massenbilanz (5.5) in Kombination mit den kinematischen Bedingungen (5.38), (5.48) und (5.70) f¨ ur die freie Oberfl¨ ache bzw. die Basis erhalten werden. Integration der Massenbilanz (5.5) von der Basis bis zur freien Oberfl¨ ache liefert unter Verwendung der Leibniz-Regel vz,s −vz,b =−∂ ∂x Zh bvxdz0+vx,s ∂h ∂x −vx,b ∂b ∂x −∂ ∂y Zh bvydz0+vy,s ∂h ∂y −vy,b ∂b ∂y.(5.110) Unter Verwendung der kinematischen Bedingungen (5.38) f¨ ur die freie Oberfl¨ ache und (5.48) bzw. (5.70) f¨ ur die Basis ergibt dies ∂H ∂t =∂(h−b) ∂t =−∂ ∂x Zh bvxdz0−∂ ∂y Zh bvydz0+a⊥ s−[ω] ˙mw b.(5.111) Diese Gleichung bilanziert die zeitliche ¨ Anderung der Eisdicke mit der horizontalen Divergenz der vertikalintegrierten Horizontalgeschwindigkeit und der AkkumulationsAblations-Funktion; letztere stellt eine klimatische Input-Gr¨ oße dar. 5.3.4 Evolution der CTS In ¨ ahnlicher Weise wie f¨ ur die freie Oberfl¨ ache kann eine Evolutionsgleichung f¨ ur die CTS hergeleitet werden. Hierzu wird die Massenbilanz (5.5) von der Basis bis zur CTS integriert, dies ergibt vz,m −vz,b =−∂ ∂x Zzm bvxdz0+vx,m ∂zm ∂x −vx,b ∂b ∂x −∂ ∂y Zzm bvydz0+vy,m ∂zm ∂y −vy,b ∂b ∂y,(5.112) 76 wobei wiederum die Leibniz-Regel verwendet wurde. Mit den kinematischen Bedingungen (5.70) f¨ ur die Basis und (5.79) f¨ ur die CTS erh¨ alt man ∂(zm−b) ∂t =−∂ ∂x Zzm bvxdz0−∂ ∂y Zzm bvydz0+a⊥ m−[ω] ˙mw b.(5.113) Die Interpretation dieses Ergebnisses ist dieselbe wie im Falle der freien Oberfl¨ ache, jedoch ist der Volumenfluß durch die CTS a⊥ mals innere Variable keine Inputgr¨ oße, so daß diese Gleichung zur Bestimmung der Evolution der CTS-Position zun¨ achst nicht geeignet ist. 5.3.5 Evolution der Eisbasis bzw. Lithosph¨ arenoberseite Hierf¨ ur ist keine weitere Rechnung vonn¨ oten, sondern (5.33) kann direkt ¨ ubernommen werden: ∂b ∂t =−1 Tr [b−(b0−ρ ρa H)].(5.114) 5.3.6 Temperatur und Wassergehalt Wie bereits erw¨ ahnt m¨ ussen die in Abschnitt 5.2 vorgestellten Gleichungen zur Berechnung des Temperaturund Wassergehaltsfeldes numerisch gel¨ ost werden. Der Vollst¨ andigkeit halber seien sie an dieser Stelle mit ihren zugeh¨ origen Randund ¨ Ubergangsbedingungen noch einmal aufgef¨ uhrt, wobei die Wassergehaltsgleichung und die CTS-¨ Ubergangsbedingung mit Hilfe des hydrostatischen Druckfeldes noch etwas vereinfacht werden k¨ onnen. Die Temperaturgleichung f¨ ur den kalten Bereich (5.13) lautet ∂θ ∂t +vx ∂θ ∂x +vy ∂θ ∂y +vz ∂θ ∂z =D c ∂ ∂z κ∂θ ∂z !+ 2α cKEA(θ0)f(σ)σ2.(5.115) Im temperierten Bereich ist die Temperatur durch den Druck eindeutig bestimmt; mit (5.28) und (5.95) erh¨ alt man θ=θM=−B(h(x, y, t)−z).(5.116) Somit gilt ∂θM ∂t +vx ∂θM ∂x +vy ∂θM ∂y +vz ∂θM ∂z =−B ∂h ∂t +vx ∂h ∂x +vy ∂h ∂y −vz!(5.117) und ∂ ∂z κ∂θM ∂z !=B∂κ ∂z =B2∂κ ∂θ ; (5.118) 77 es ergibt sich f¨ ur die Wassergehaltsgleichung im temperierten Bereich (5.30) ∂ω ∂t +vx ∂ω ∂x +vy ∂ω ∂y +vz ∂ω ∂z =Dt ∂2ω ∂z2+Dαt αB2∂κ ∂θ +cαt αB ∂h ∂t +vx ∂h ∂x +vy ∂h ∂y −vz!+ 2αtKEAt(ω)ft(σ)σ2,(5.119) wobei f¨ ur die effektive Schubspannung σgem¨ aß (5.98) σ= (h−z)v u u t ∂h ∂x!2 + ∂h ∂y !2 (5.120) gilt. Die Temperaturgleichung f¨ ur die Lithosph¨ are (5.32) lautet mit (5.34) – (5.36) ∂θ ∂t +∂b ∂t ∂θ ∂z =Drκr cr ∂2θ ∂z2.(5.121) An der kalten freien Oberfl¨ ache ist die Temperatur vorgeschrieben (Gl. (5.47), Randbedingung vom Dirichlet-Typ): θ=θs(x, y, t); (5.122) im unwahrscheinlichen Fall einer temperierten freien Oberfl¨ ache wird statt dessen der Wassergehalt ωoder dessen Normalableitung als Randbedingung verwendet. F¨ ur eine kalte Eisbasis gelten (5.67) und (5.68): κ∂θ ∂z −[κr] [κ]κr ∂θ+ ∂z =−α D[(vsl)xσxz + (vsl)yσyz],(5.123) θ=θ+,(5.124) und beim Vorliegen einer temperierten Basis kommen (5.72), (5.74), (5.75) zur Anwendung: Dt ∂ω ∂z = (1 −[ω]ω) ˙mw b−Pw b,(5.125) θ=θ+=θM,(5.126) Pw b=Dαt ακ∂θ ∂z −Dαt α [κr] [κ]κr ∂θ+ ∂z +αt((vsl)xσxz + (vsl)yσyz).(5.127) Gleichung (5.125) wird allerdings nur ben¨ otigt, wenn der temperierten Basis eine temperierte Eisschicht nichtverschwindender Dicke ¨ uberlagert ist. In diesem Fall kann (5.127) noch vereinfacht werden, da dann wegen (5.116) ∂θ/∂z =Bgilt. F¨ ur den Fall einer temperierten Basis, die von kaltem Eis ¨ uberlagert ist, ist dies nat¨ urlich nicht m¨ oglich. 78 An der Lithosph¨ arenunterseite gilt (5.77): κr ∂θ− ∂z =−NrQ⊥ geoth.(5.128) Es verbleiben noch die ¨ Ubergangsbedingungen an der CTS (5.83), (5.89), (5.91) und (5.93): θ+=θ−=θM.(5.129) i) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)>0 (Schmelzbedingung): ω−= 0,(5.130) ∂θ+ ∂z =B.(5.131) ii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)<0 (Gefrierbedingung): Dκ ∂θ+ ∂z −B!−α αtDt ∂ω− ∂z =α αt ω−a⊥ m.(5.132) Hierf¨ ur wurde ∂θ− M/∂z =Bverwendet. iii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z) = 0 (Grenzfall): ω−≥0 (unbestimmt),(5.133) ∂θ+ ∂z =B.(5.134) Abschließend sei auf einige Beschr¨ ankungen der SIA-N¨ aherung hingewiesen. Zun¨ achst schließt diese Skalierung geschlossene CTS-Linien aus, da sie nur flache CTS-Linien erlaubt. Einschl¨ usse von temperiertem Eis in einer Umgebung von kaltem Eis sind also nicht erfaßt, sondern nur solche temperierten Bereiche, die bis zur Eisbasis reichen. Das ist aber der weitaus wichtigste Fall. Weiterhin ergibt die SIA Singularit¨ aten an Eisr¨ andern und an Eisdomen (Maxima der Eish¨ ohe ¨ uber Normal Null). Hutter [31] weist darauf hin, daß unter der Voraussetzung einer beschr¨ ankten basalen Gleitfunktion C(t)die Eish¨ ohe an seitlichen Begrenzungen in Form einer Wurzelfunktion gegen Null strebt und somit am Rand eine vertikale Tangente aufweist, was im Widerspruch zur Flachheitsannahme steht. F¨ ur Eisdome erh¨ alt man bei Verwendung einer pseudoplastischen Kriechfunktion f(t)(σ) = σn−1mit n > 1, wie sie ¨ ublicherweise angesetzt wird (vgl. §5.5), eine unendliche Kr¨ ummung der Eisoberfl¨ ache, was mit 79 der gegen unendlich strebenden Viskosit¨ at einer solchen Kriechfunktion bei kleinen Spannungen bzw. Verzerrungsgeschwindigkeiten zusammenh¨ angt (Hutter, Yakowitz & Szidarovszky [33]). Diese Ergebnisse sind f¨ ur den Fall eines reinen Kalteismodells hergeleitet, jedoch auf die SIA f¨ ur polytherme Eisschilde ¨ ubertragbar. Weiterhin ist bei Eisdomen die SIA-Annahme, daß die auftretenden Vertikalgeschwindigkeiten von der Ordnung Aspektverh¨ altnis ×Horizontalgeschwindigkeit sind, nicht mehr realistisch, da dort die SIA-Approximation ein rein vertikales Geschwindigkeitsfeld ergibt; vgl. Gleichungen (5.107), (5.108) und (5.109). Eine strenge Behandlung m¨ ußte daher in der N¨ ahe von Eisdomen und seitlichen R¨ andern die SIA aufgeben. Es kann jedoch davon ausgegangen werden, daß eine numerische L¨ osung des Modells mit diskreten Gitterpunkten diese lokalen Singularit¨ aten ausreichend verschmiert, so daß sich der verbleibende Fehler in vertretbaren Grenzen h¨ alt. 5.4 Zusammenstellung in dimensionsbehafteter Form Der ¨ Ubersichtlichkeit und besseren Anschaulichkeit halber sollen an dieser Stelle die in Abschnitt 5.3 hergeleiteten teilintegrierten polythermen SIA-Gleichungen in dimensionsbehafteter Form aufgef¨ uhrt werden. Diese lassen sich leicht aus den dimensionslosen Gleichungen gewinnen, indem alle typischen Werte [·] auf 1 gesetzt werden (vgl. Abschnitt 5.1). Spannungen: p=ρg(h−z),(5.135) σxz =−ρg(h−z)∂h ∂x,(5.136) σyz =−ρg(h−z)∂h ∂y ,(5.137) σ=ρg(h−z)v u u t ∂h ∂x!2 + ∂h ∂y !2 .(5.138) Geschwindigkeit: vx=−ρg(h−b)C(t) ∂h ∂x −2ρg∂h ∂x Zz bEA(t)(·)f(t)(σ)(h−z0)dz0,(5.139) vy=−ρg(h−b)C(t) ∂h ∂y −2ρg∂h ∂y Zz bEA(t)(·)f(t)(σ)(h−z0)dz0,(5.140) vz=−Zz b ∂vx ∂x +∂vy ∂y !dz0+∂b ∂t +vx,b ∂b ∂x +vy,b ∂b ∂y −˙mw b ρ.(5.141) 80 Entwicklung der freien Oberfl¨ ache: ∂H ∂t =∂(h−b) ∂t =−∂qx ∂x −∂qy ∂y +a⊥ s−˙mw b ρ,(5.142) mit dem neu eingef¨ uhrten Massenfluß q, definiert als (qx, qy) := Zh b(vx, vy)dz0.(5.143) Entwicklung der CTS: ∂(zm−b) ∂t =−∂ ∂x Zzm bvxdz0−∂ ∂y Zzm bvydz0+a⊥ m−˙mw b ρ.(5.144) Entwicklung der Eisbasis bzw. Lithosph¨ arenoberseite: ∂b ∂t =−1 τV [b−(b0−ρ ρa H)].(5.145) Temperatur und Wassergehalt: Temperaturgleichung, kalter Bereich: ∂T ∂t +vx ∂T ∂x +vy ∂T ∂y +vz ∂T ∂z =1 ρc ∂ ∂z κ∂T ∂z !+2 ρcEA(T0)f(σ)σ2.(5.146) Temperatur, temperierter Bereich: T=TM=T0−β(h−z).(5.147) Wassergehaltsgleichung, temperierter Bereich: ∂ω ∂t +vx ∂ω ∂x +vy ∂ω ∂y +vz ∂ω ∂z =ν ρ ∂2ω ∂z2+β2 ρL ∂κ ∂T +cβ L ∂h ∂t +vx ∂h ∂x +vy ∂h ∂y −vz!+2 ρLEAt(ω)ft(σ)σ2−1 ρD(ω).(5.148) Hier wurde ein zus¨ atzlicher Term −D(ω)/ρ eingef¨ uhrt; D(ω) ist die WasserdrainageFunktion, die einer negativen volumenm¨ aßigen Wasserproduktion entspricht und einen einfachen ad-hoc-Ansatz der Beschreibung der Wasserdrainage aus den temperierten Eisbereichen in den Boden hinein darstellt. Dies ist notwendig, da, wie sich bei vorl¨ aufigen numerischen Rechnungen ohne Drainage gezeigt hat, ansonsten unrealistische und sogar unphysikalische (ω > 100%) Werte f¨ ur den Wassergehalt resultieren k¨ onnen. Es wird in Kauf genommen, daß dieses “Wegbeamen” von Wasser aus dem Eisinnern eine Verletzung der lokalen Massenbilanz darstellt. Realistischere Drainage-Mechanismen, die auf dem Zusammenspiel von Schwerkraft und Wechselwirkungskr¨ aften zwischen Wasser und Eis beruhen, w¨ urden ein noch komplizierteres 81 Modell mit zwei getrennten Impulsbilanzen f¨ ur Wasser und Eis erforderlich machen, was hier nicht verfolgt werden soll (vgl. Bauer [3], Wu [70, 71]). Temperaturgleichung, Lithosph¨ are: ∂T ∂t +∂b ∂t ∂T ∂z =κr ρrcr ∂2T ∂z2.(5.149) Kalte freie Oberfl¨ ache: T=Ts(x, y, t).(5.150) Kalte Eisbasis: κ∂T ∂z −κr ∂T+ ∂z =−(vsl)xσxz −(vsl)yσyz,(5.151) T=T+.(5.152) Temperierte Eisbasis (keine temperierte Schicht dar¨ uber): T=T+=TM,(5.153) Pw b=1 L(κ∂T ∂z −κr ∂T+ ∂z + ((vsl)xσxz + (vsl)yσyz)).(5.154) Temperierte Eisbasis (von temperierter Schicht ¨ uberlagert): ν∂ω ∂z = (1 −ω) ˙mw b−Pw b,(5.155) T=T+=TM,(5.156) Pw b=1 L(κβ −κr ∂T+ ∂z + ((vsl)xσxz + (vsl)yσyz)).(5.157) Es sei an dieser Stelle noch einmal bemerkt, daß im Falle vernachl¨ assigbarer Wasserdiffusivit¨ at νder im allgemeinen vorzugebende Wasser-Massenstrom in den Boden ˙mw bvia (3.48) durch die basale Schmelzrate Pw bausgedr¨ uckt werden kann. Da der Wassergehalt im temperierten Eis als klein angenommen wird (vgl. auch (5.173)), gilt folglich ˙mw b≈ Pw b. Lithosph¨ arenunterseite: κr ∂T− ∂z =−Q⊥ geoth.(5.158) ¨ Ubergangsbedingungen, CTS: T+=T−=TM.(5.159) i) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)>0 (Schmelzbedingung): ω−= 0,(5.160) 82 ∂T+ ∂z =β. (5.161) ii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)<0 (Gefrierbedingung): κ ∂T+ ∂z −β!−Lν ∂ω− ∂z =Lρω−a⊥ m.(5.162) iii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z) = 0 (Grenzfall): ω−≥0 (unbestimmt),(5.163) ∂T+ ∂z =β. (5.164) Es sei daran erinnert (vgl. §3.2.5), daß die im allgemeinen Fall aufwendige Unterscheidung nach Schmelzbzw. Gefrierbedingungen im Falle vernachl¨ assigbarer Wasserdiffusion einfach nach dem Vorzeichen des Volumenflusses durch die CTS, a⊥ m, erfolgt. 5.5 Spezifizierung physikalischer Gr¨ oßen Die im vorigen Abschnitt zusammengestellten Gleichungen des polythermen SIAEismodells enthalten eine Reihe von noch unbestimmten physikalischen Gr¨ oßen. Diese sollen hier spezifiziert werden. Im Rahmen dieser Arbeit werden im einzelnen verwendet: Kriechfunktion f¨ ur kaltes und temperiertes Eis: Glen’sches Fließgesetz (vgl. Glen [19], Nye [57], Hooke [26], Paterson [62]) f(σ) = ft(σ) = σn−1mit n= 3.(5.165) Rate-Faktor f¨ ur kaltes Eis: Arrhenius-Gesetz (genaue Form nach Paterson [62]) A(T0) = A0e−Q/R(T0+T0),(5.166) mit der Aktivierungsenergie Q Q= 60 kJ mol−1(T0<−10◦C), Q= 139 kJ mol−1(T0≥ −10◦C),(5.167) der universellen Gaskonstanten R, und als verbindendem Wert f¨ ur die beiden Temperaturbereiche A(T0=−10◦C) = 5.2·10−25 s−1Pa−3.(5.168) 83 − ∂b ∂x +ζ0 t ∂Ht ∂x !∂vx,t ∂ζ0 t− ∂b ∂y +ζ0 t ∂Ht ∂y !∂vy,t ∂ζ0 t)dζ0 t +∂b ∂t +vx,b ∂b ∂x +vy,b ∂b ∂y −˙mw b ρ.(6.23) Geschwindigkeit, kalter Bereich: vx,c =vx(zm)−2(ρg)3kgrad hk2∂h ∂x H4 ca ea−1× Zζc 0EA(T0 c)(1 −E(ζ0 c))3eaζ0 cdζ0 c,(6.24) vy,c =vy(zm)−2(ρg)3kgrad hk2∂h ∂y H4 ca ea−1× Zζc 0EA(T0 c)(1 −E(ζ0 c))3eaζ0 cdζ0 c,(6.25) vz,c =vz(zm)−Zζc 0(Hcaeaζ0 c ea−1 ∂vx,c ∂ξc +∂vy,c ∂ηc! − ∂zm ∂x +E(ζ0 c)∂Hc ∂x !∂vx,c ∂ζ0 c− ∂zm ∂y +E(ζ0 c)∂Hc ∂y !∂vy,c ∂ζ0 c)dζ0 c.(6.26) Entwicklung der freien Oberfl¨ ache: ∂H ∂t =∂(h−b) ∂t =−∂qx ∂x −∂qy ∂y +a⊥ s−˙mw b ρ,(6.27) mit (qx, qy) = HtZ1 0(vx,t, vy,t)dζ0 t+Hca ea−1Z1 0(vx,c, vy,c)eaζ0 cdζ0 c.(6.28) Entwicklung der CTS: ∂Ht ∂t =−∂ ∂x HtZ1 0vx,t dζ0 t−∂ ∂y HtZ1 0vy,t dζ0 t+a⊥ m−˙mw b ρ.(6.29) Entwicklung der Eisbasis bzw. Lithosph¨ arenoberseite: ∂b ∂t =−1 τV [b−(b0−ρ ρa H)].(6.30) Temperatur und Wassergehalt: Temperaturgleichung, kalter Bereich: ∂Tc ∂τc +vx,c ∂Tc ∂ξc +vy,c ∂Tc ∂ηc +∂Tc ∂ζc(vz,c Hcaeaζc(ea−1) −1 Hcaeaζc (ea−1)∂zm ∂t + (eaζc−1)∂Hc ∂t ! −vx,c Hcaeaζc (ea−1)∂zm ∂x + (eaζc−1)∂Hc ∂x ! 90 −vy,c Hcaeaζc (ea−1)∂zm ∂y + (eaζc−1)∂Hc ∂y !) =ea−1 ρc(Tc)Hcaeaζc ∂ ∂ζc κ(Tc)ea−1 Hcaeaζc ∂Tc ∂ζc! +2 g c(Tc)(ρg)3EA(T0 c)H4 c(1 −E(ζc))4kgrad hk4.(6.31) Schmelztemperatur, kalter Bereich: TM,c =T0−βHc(1 −E(ζc)).(6.32) Temperatur, temperierter Bereich: Tt=TM,t =T0−β(Hc+Ht(1 −ζt)).(6.33) Wassergehaltsgleichung, temperierter Bereich: ∂ωt ∂τt +vx,t ∂ωt ∂ξt +vy,t ∂ωt ∂ηt +∂ωt ∂ζt(vz,t Ht−1 Ht ∂b ∂t +ζt ∂Ht ∂t !−vx,t Ht ∂b ∂x +ζt ∂Ht ∂x !−vy,t Ht ∂b ∂y +ζt ∂Ht ∂y !) =ν ρH2 t ∂2ωt ∂ζ2 t + 2 g L(ρg)3EAt(ωt)(Hc+Ht(1 −ζt))4kgrad hk4 +β2 ρL ∂κ ∂T (TM,t) + c(TM,t)β L ∂h ∂t +vx,t ∂h ∂x +vy,t ∂h ∂y −vz,t!−D(ωt) ρ.(6.34) Temperaturgleichung, Lithosph¨ are: ∂Tr ∂τr +∂Tr ∂ζr(1 Hr ∂b ∂t −1 Hr ∂br ∂t +ζr ∂Hr ∂t !)=κr ρrcrH2 r ∂2Tr ∂ζ2 r .(6.35) Kalte freie Oberfl¨ ache: Tc=Ts(x, y, t).(6.36) Kalte Eisbasis: Tc=Tr,(6.37) κ(Tc)ea−1 Hca ∂Tc ∂ζc−κr Hr ∂Tr ∂ζr = (vsl)xρgHc ∂h ∂x + (vsl)yρgHc ∂h ∂y .(6.38) Temperierte Eisbasis (keine temperierte Schicht dar¨ uber): Tc=Tr=TM,c,(6.39) Pw b=1 L(κ(Tc)ea−1 Hca ∂Tc ∂ζc−κr Hr ∂Tr ∂ζr −(vsl)xρgHc ∂h ∂x −(vsl)yρgHc ∂h ∂y ).(6.40) 91 Temperierte Eisbasis (von temperierter Schicht ¨ uberlagert): Tt=Tr=TM,t,(6.41) ν Ht ∂ωt ∂ζt = (1 −ωt) ˙mw b−Pw b,(6.42) Pw b=1 L(κ(Tt)β−κr Hr ∂Tr ∂ζr−(vsl)xρgH ∂h ∂x −(vsl)yρgH ∂h ∂y ).(6.43) Lithosph¨ arenunterseite: κr Hr ∂Tr ∂ζr =−Q⊥ geoth.(6.44) ¨ Ubergangsbedingungen, CTS: Tc=Tt=TM,t.(6.45) i) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)>0 (Schmelzbedingung): ωt= 0,(6.46) ea−1 Hca ∂Tc ∂ζc =β. (6.47) ii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z)<0 (Gefrierbedingung): κ(Tc) ea−1 Hca ∂Tc ∂ζc−β!−Lν Ht ∂ωt ∂ζt =Lρωta⊥ m.(6.48) iii) −(wx−vw,x)∂zm ∂x −(wy−vw,y)∂zm ∂y + (wz−vw,z) = 0 (Grenzfall): ωt≥0 (unbestimmt),(6.49) ea−1 Hca ∂Tc ∂ζc =β. (6.50) Dieser Fall ist bei numerischen Berechnungen nicht von Bedeutung, da exakte Gleichheit beim Rechnen mit Gleitkommazahlen nicht auftritt. Er wird daher im folgenden nicht mehr betrachtet. 6.3 Das numerische Gitter Zur numerischen L¨ osung der auf σ-Koordinaten transformierten Modellgleichungen des vorigen Abschnitts mit einem Finite-Differenzen-Verfahren wird zun¨ achst ein numerisches Gitter eingef¨ uhrt. Dieses wird im Raum der σ-Koordinaten ¨ aquidistant 92 gew¨ ahlt; die Diskretisierung des Raumes und der Zeit ist daher x(i) = ξc(i) = ξt(i) = ξr(i) = x0+i∆x(i= 0 . . . imax), y(j) = ηc(j) = ηt(j) = ηr(j) = y0+j∆y(j= 0 . . . jmax), ζc(kc) = kc/kc,max (kc= 0 . . . kc,max), ζt(kt) = kt/kt,max (kt= 0 . . . kt,max), ζr(kr) = kr/kr,max (kr= 0 . . . kr,max), t(n) = τc(n) = τt(n) = τr(n) = t0+n∆t(n= 0 . . . nmax), t(˜n) = τc(˜n) = τt(˜n) = τr(˜n) = t0+ ˜nf ∆t(˜n= 0 . . . ˜nmax). (6.51) Hierbei bedeuten ∆xund ∆ydie Aufl¨ osungen in der Horizontalen; die Vertikalaufl¨ osungen f¨ ur die drei Bereiche kaltes Eis, temperiertes Eis und Lithosph¨ are sind ∆ζc,t,r = 1/kc,t,r,max. Zwei verschiedene Zeitschrittweiten, n¨ amlich ∆tund f ∆t, treten auf, wobei f ∆tein ganzzahliges Vielfaches von ∆tist. Diese unterschiedlichen Zeitschrittweiten werden verwendet, weil die Berechnung der Temperatur, des Wassergehaltes, des Eisalters und die Positionierung der CTS mit dem gr¨ oßeren f ∆tdurchgef¨ uhrt werden, alle ¨ ubrigen Berechnungen jedoch mit dem kleineren ∆t. Dies dient dazu, Rechenzeit bei den Simulationen einzusparen, welche bei 3-d-Str¨ omungsproblemen stets ein kritischer Punkt ist. Abbildung 6.2: Vertikalschnitt durch den Kalteisbereich des Rechengebietes. Die Geschwindigkeiten sind zwischen den Gitterpunkten definiert (Staggered Grid). 93 Zwecks besserer numerischer Stabilit¨ at und Uniformit¨ at in der Approximation wird ein Staggered Grid (gestaffeltes Gitter) verwendet, wobei die drei Komponenten der Geschwindigkeit vim kalten bzw. temperierten Eisbereich, die zwei Komponenten des Massenflusses qsowie die jeweils zwei Komponenten der Gradienten von h,zm, b,br,Hc,Htund Hrnicht auf, sondern zwischen benachbarten Gitterpunkten der betreffenden Richtung berechnet werden (vgl. auch Abbildung 6.2): vx,c(i±1 2, j, kc, n), vy,c(i, j ±1 2, kc, n), vz,c(i, j, kc±1 2, n), vx,t(i±1 2, j, kt, n), vy,t(i, j ±1 2, kt, n), vz,t(i, j, kt±1 2, n), qx(i±1 2, j, n), qy(i, j ±1 2, n), ∂h, zm, b, br, Hc, Ht, Hr ∂x (i±1 2, j, n),∂h, zm, b, br, Hc, Ht, Hr ∂y (i, j ±1 2, n). Alle ¨ ubrigen Gr¨ oßen sind auf den Gitterpunkten selbst definiert, z. B. die Temperatur im kalten Eis Tc(i, j, kc,˜n), oder der Wassergehalt im temperierten Bereich ωt(i, j, kt,˜n). Zuweilen werden jedoch bei der numerischen L¨ osung Gr¨ oßen an Stellen ben¨ otigt, wo sie eigentlich nicht definiert sind; in diesem Fall wird lineare Interpolation angewendet, um sie dorthin zu transportieren. Dies wird im folgenden durch einen Querstrich ¨ uber der betreffenden Gr¨ oße kenntlich gemacht. Beispielsweise ist ¯vx,c(i, j, kc, n) := 1 2vx,c(i−1 2, j, kc, n) + vx,c(i+1 2, j, kc, n), ∂h ∂x(i, j, n) := 1 2 ∂h ∂x(i−1 2, j, n) + ∂h ∂x(i+1 2, j, n)!, ¯ Tc(i, j, kc+1 2,˜n) := 1 2Tc(i, j, kc,˜n) + Tc(i, j, kc+ 1,˜n).(6.52) Die Notwendigkeit solcher Interpolationen, welche zun¨ achst sehr umst¨ andlich anmuten, ist gerade der große Vorteil des Staggered Grid, denn durch diese Mittelungen werden sich entwickelnde numerische Instabilit¨ aten recht wirkungsvoll ged¨ ampft. 6.4 Diskretisierung der Modellgleichungen Im folgenden wird das den numerischen Berechnungen dieser Arbeit zugrundeliegende Diskretisierungsschema f¨ ur die Modellgleichungen angegeben. Die Reihenfolge ist dabei identisch mit der im Computerprogramm tats¨ achlich verwendeten Abarbeitungsfolge im Rahmen einer zeitlichen Iteration. Das Problem der Bestimmung der CTS-Position wird zun¨ achst noch ausgeklammert; dies wird im folgenden Abschnitt ausf¨ uhrlich besprochen. Im allgemeinen muß zwischen drei verschiedenen F¨ allen unterschieden werden, n¨ amlich einer kalten Eisbasis, einer temperierten Eisbasis ohne ¨ uberlagerte tempe94 rierte Eisschicht, und einer temperierten Eisbasis mit ¨ uberlagerter temperierter Eisschicht. Hat die Eiss¨ aule ¨ uber dem horizontalen Gitterpunkt (i, j) eine kalte Basis, soll dies als Fall CB bezeichnet werden, bei einer temperierten Basis ohne ¨ uberlagerte temperierte Schicht als Fall TB, und bei einer temperierten Basis mit temperierter Schicht als Fall TL. Zur Vereinfachung der Notation werden die Zeitindices nbzw. ˜nnicht mehr explizit aufgef¨ uhrt. Statt dessen werden Gr¨ oßen, die bei einem zeitlichen Iterationsschritt von z. B. n0nach n0+ 1 f¨ ur den neuen Zeitpunkt n0+ 1 gelten sollen, mit einem hochgestellten Stern gekennzeichnet. h(i, j) ist in diesem Fall also gleichbedeutend mit h(i, j, n0), h?(i, j) mit h(i, j, n0+ 1); entsprechendes gilt f¨ ur alle ¨ ubrigen Gr¨ oßen. Temperatur und Wassergehalt Die Gleichungen (6.31), (6.34) und (6.35) enthalten jeweils die zeitliche Ableitung der Temperatur bzw. des Wassergehaltes, sind also prognostischer Natur. Wie bereits weiter oben erw¨ ahnt, wird das Problem der Bestimmung der CTS-Position hier noch nicht besprochen. Die Diskretisierungen werden daher unter der Annahme angegeben, daß bereits bekannt sei, ob Fall CB, Fall TB oder Fall TL vorliegt, und daß bei letzterem die CTS-Position bestimmt sei. Die Diskretisierung wird derart durchgef¨ uhrt, daß Ableitungen in vertikaler Richtung implizit gemacht werden, also zum jeweils neuen Zeitpunkt genommen werden; andere Ableitungen werden dagegen explizit diskretisiert. Dies erm¨ oglicht bereits deutlich gr¨ oßere Zeitschrittweiten verglichen mit einem voll expliziten Verfahren, da bei tats¨ achlichen Simulationen in vertikaler Richtung wesentlich kleinere Gitterweiten auftreten als in horizontaler Richtung. Andererseits sind durch die Beschr¨ ankung auf implizite Vertikalableitungen lediglich lineare Gleichungssysteme mit tridiagonaler Matrix (f¨ ur jede Vertikals¨ aule des Rechengebietes getrennt) zu l¨ osen, was schnell und unproblematisch durchgef¨ uhrt werden kann. F¨ ur die horizontalen Advektionsterme in (6.31), (6.34) wird zur numerischen Stabilisierung ein Upstream- (auch: Upwind-) Diskretisierungsschema verwendet, wobei Horizontalgeschwindigkeiten und horizontale Temperaturbzw. Wassergehaltsableitungen um einen halben Gitterpunkt in xbzw. y-Richtung versetzt genommen werden. Diese Verschiebung erfolgt stets stromaufw¨ arts, also entgegen der Richtung von vxbzw. vy; sie kann auch als numerische Diffusion ausgedr¨ uckt werden. Fall CB: Lithosph¨ arenunterseite: κr Hr·T? r(i, j, 1) −T? r(i, j, 0) ∆ζr =−Q⊥ geoth.(6.53) 95 Temperaturgleichung, Lithosph¨ are: T? r(i, j, kr)−Tr(i, j, kr) f ∆t =κr ρrcrH2 r·T? r(i, j, kr+ 1) −2T? r(i, j, kr) + T? r(i, j, kr−1) ∆ζ2 r (kr= 1 . . . kr,max −1).(6.54) Kalte Eisbasis: T? c(i, j, 0) = T? r(i, j, kr,max),(6.55) κ(Tc(i, j, 0)) (ea−1) aHc(i, j)·T? c(i, j, 1) −T? c(i, j, 0) ∆ζc −κr Hr·T? r(i, j, kr,max)−T? r(i, j, kr,max −1) ∆ζr =ρg ¯vx,c(i, j, 0) Hc(i, j)∂h ∂x(i, j) +ρg ¯vy,c(i, j, 0) Hc(i, j)∂h ∂y (i, j).(6.56) Temperaturgleichung, kalter Bereich: T? c(i, j, kc)−Tc(i, j, kc) f ∆t +1 2vx,c(i+1 2, j, kc)−|vx,c(i+1 2, j, kc)|·Tc(i+ 1, j, kc)−Tc(i, j, kc) ∆x +1 2vx,c(i−1 2, j, kc) + |vx,c(i−1 2, j, kc)|·Tc(i, j, kc)−Tc(i−1, j, kc) ∆x +1 2vy,c(i, j +1 2, kc)−|vy,c(i, j +1 2, kc)|·Tc(i, j + 1, kc)−Tc(i, j, kc) ∆y +1 2vy,c(i, j −1 2, kc) + |vy,c(i, j −1 2, kc)|·Tc(i, j, kc)−Tc(i, j −1, kc) ∆y +T? c(i, j, kc+ 1) −T? c(i, j, kc−1) 2∆ζc× (¯vz,c(i, j, kc) aHc(i, j)eaζc(kc)(ea−1) −1 aHc(i, j)eaζc(kc) (ea−1) ∂zm ∂t (i, j)+(eaζc(kc)−1) ∂Hc ∂t (i, j)! −¯vx,c(i, j, kc) aHc(i, j)eaζc(kc) (ea−1) ∂zm ∂x (i, j)+(eaζc(kc)−1) ∂Hc ∂x (i, j)! −¯vy,c(i, j, kc) aHc(i, j)eaζc(kc) (ea−1) ∂zm ∂y (i, j)+(eaζc(kc)−1) ∂Hc ∂y (i, j)!) =ea−1 aρc(Tc(i, j, kc)) Hc(i, j)eaζc(kc)∆ζc× 96      κ(¯ Tc(i, j, kc+1 2)) (ea−1) a¯ Hc(i, j, kc+1 2)ea¯ ζc(kc+1 2)·T? c(i, j, kc+ 1) −T? c(i, j, kc) ∆ζc −κ(¯ Tc(i, j, kc−1 2)) (ea−1) a¯ Hc(i, j, kc−1 2)ea¯ ζc(kc−1 2)·T? c(i, j, kc)−T? c(i, j, kc−1) ∆ζc     + 2 g c(Tc(i, j, kc))(ρg)3EA(T0 c(i, j, kc)) H4 c(i, j)1−E(ζc(kc))4×   ∂h ∂x!2 (i, j) + ∂h ∂y !2 (i, j)  2 (kc= 1 . . . kc,max −1).(6.57) Kalte freie Oberfl¨ ache: T? c(i, j, kc,max) = Ts(x(i), y(j), t?).(6.58) Lithosph¨ arentemperatur und Kalteistemperatur stellen ein gekoppeltes Problem dar und m¨ ussen simultan berechnet werden. Durch Sortieren der Terme nach den Unbekannten T? r(i, j, kr) und T? c(i, j, kc) erh¨ alt man tridiagonale Gleichungssysteme des Ranges kr,max +kc,max +2, die unabh¨ angig voneinander f¨ ur alle S¨ aulen (i, j) mit kalter Basis zu l¨ osen sind. Gleichung (6.54) enth¨ alt, verglichen mit der urspr¨ unglichen Gleichung (6.35), keine erste Ableitung der Temperatur nach ζrmehr. Dies kommt daher, daß f¨ ur die in dieser Arbeit durchgef¨ uhrten Simulationen stets mit konstanter Lithosph¨ arendicke Hr(vgl. §5.5) gerechnet wird, so daß die eingeklammerte Differenz in (6.35) verschwindet. Der allgemeine Fall ist jedoch problemlos zu implementieren. Fall TB: Dies entspricht weitgehend dem Fall CB; lediglich (6.55), (6.56) sind zu ersetzen durch das Pendant f¨ ur die temperierte Eisbasis: T? c(i, j, 0) = TM,c(i, j, 0),(6.59) T? r(i, j, kr,max) = TM,c(i, j, 0).(6.60) Als Folge hiervon entkoppeln die Berechnungen von Lithosph¨ arentemperatur und Kalteistemperatur; man erh¨ alt getrennte tridiagonale Gleichungssysteme f¨ ur T? r(i, j, kr) (vom Rang kr,max + 1) und T? c(i, j, kc) (vom Rang kc,max + 1). Fall TL: Auch hier sind Lithosph¨ are und Eisbereich entkoppelt. Die Berechnung der Lithosph¨ arentemperatur erfolgt analog zum Fall TB nach (6.53), (6.54) und (6.60) und ergibt wiederum tridiagonale Gleichungssysteme f¨ ur T? r(i, j, kr) vom Rang kr,max + 1. 97 F¨ ur die Berechnung des Wassergehaltes im temperierten Bereich wird zun¨ achst ausgenutzt, daß keine physikalische Wasserdiffusion ber¨ ucksichtigt werden soll (§5.5, Tabelle 5.1). Daher verschwindet (6.42) identisch, der basale Wassergehalt ist also keiner Randbedingung mehr unterworfen. Weiterhin kann an der CTS die Unterscheidung nach Schmelzbzw. Gefrierbedingungen einfach nach dem Vorzeichen des Volumenstromes durch die CTS a⊥ mvorgenommen werden, und der Diffusionsterm in (6.48) verschwindet ebenfalls. Zwecks numerischer Stabilit¨ at wird jedoch der vertikale Diffusionsterm in der Wassergehaltsgleichung (6.34) beibehalten und als numerische Diffusion interpretiert, wobei die kleine Diffusivit¨ at ν= 10−6kg m−1s−1verwendet wird. Dies f¨ uhrt aber wegen der parabolischen Natur von Diffusionsgleichungen dazu, daß f¨ ur den basalen Wassergehalt eine Randbedingung n¨ otig ist, welche die Physik nicht zur Verf¨ ugung stellt. Um das Ergebnis m¨ oglichst wenig zu verf¨ alschen, wird eine Neumannsche Bedingung angenommen, n¨ amlich das Verschwinden der Vertikalableitung des Wassergehaltes an der Basis. Ein ¨ ahnliches Problem tritt an der CTS als Obergrenze des temperierten Eisbereiches auf. Im Falle von Schmelzbedingungen gilt die Dirichlet-Bedingung (6.46); f¨ ur Gefrierbedingungen ist jedoch keine Randbedingung vorhanden (Gleichung (6.48) wird zur Berechnung des Temperaturfeldes im Kalteisbereich verwendet, siehe weiter unten, und kann daher nicht auch f¨ ur die Berechnung des Wassergehaltes ausgenutzt werden), so daß wiederum das Verschwinden der Vertikalableitung des Wassergehaltes als numerische Randbedingung angesetzt wird4. Temperierte Eisbasis: numerische Randbedingung ω? t(i, j, 1) −ω? t(i, j, 0) ∆ζt = 0.(6.61) Wassergehaltsgleichung, temperierter Bereich: ω? t(i, j, kt)−ωt(i, j, kt) f ∆t +1 2vx,t(i+1 2, j, kt)−|vx,t(i+1 2, j, kt)|·ωt(i+ 1, j, kt)−ωt(i, j, kt) ∆x +1 2vx,t(i−1 2, j, kt) + |vx,t(i−1 2, j, kt)|·ωt(i, j, kt)−ωt(i−1, j, kt) ∆x +1 2vy,t(i, j +1 2, kt)−|vy,t(i, j +1 2, kt)|·ωt(i, j + 1, kt)−ωt(i, j, kt) ∆y 4Nat¨ urlich birgt die Einf¨ uhrung von solchen k¨ unstlichen Randbedingungen eine gewisse Problematik, da sie die Resultate beeinflussen, auch wenn die Verf¨ alschung durch eine hinreichend kleine Diffusivit¨ at und Wahl einer Neumannschen statt Dirichletschen Bedingung begrenzt werden kann. Leider besteht bei numerischen Simulationen h¨ aufig die Notwendigkeit, derartige Kompromisse einzugehen, um ¨ uberhaupt stabile Modell¨ aufe zu erm¨ oglichen. 98 +1 2vy,t(i, j −1 2, kt) + |vy,t(i, j −1 2, kt)|·ωt(i, j, kt)−ωt(i, j −1, kt) ∆y +ω? t(i, j, kt+ 1) −ω? t(i, j, kt−1) 2∆ζt× (¯vz,t(i, j, kt) H(?) t(i, j)−1 H(?) t(i, j) ∂b ∂t(i, j) + ζt(kt)∂Ht ∂t (?) (i, j)! −¯vx,t(i, j, kt) H(?) t(i, j) ∂b ∂x(i, j) + ζt(kt)∂Ht ∂x (i, j)! −¯vy,t(i, j, kt) H(?) t(i, j) ∂b ∂y(i, j) + ζt(kt)∂Ht ∂y (i, j)!) =ν ρH(?) t2(i, j)·ω? t(i, j, kt+ 1) −2ω? t(i, j, kt) + ω? t(i, j, kt−1) ∆ζ2 t + 2 g L(ρg)3EAt(ωt(i, j, kt)) H(?) c(i, j) + H(?) t(i, j) (1 −ζt(kt))4×   ∂h ∂x!2 (i, j) + ∂h ∂y !2 (i, j)  2 +β2 ρL ∂κ ∂T +βc(TM,t(i, j, kt)) L× ∂h ∂t (i, j) + ¯vx,t(i, j, kt)∂h ∂x(i, j) +¯vy,t(i, j, kt)∂h ∂y (i, j)−¯vz,t(i, j, kt)!(kt= 1 . . . kt,max −1).(6.62) CTS, a⊥ m(?)(i, j)>0: physikalische Randbedingung ω? t(i, j, kt,max) = 0.(6.63) CTS, a⊥ m(?)(i, j)<0: numerische Randbedingung ω? t(i, j, kt,max)−ω? t(i, j, kt,max −1) ∆ζt = 0.(6.64) In dieser Formulierung ist die Berechnung des Wassergehaltes im temperierten Eisbereich entkoppelt von Temperaturberechnungen in der Lithosph¨ are bzw. im Kalteisbereich; das Problem besteht aus tridiagonalen Gleichungssystemen f¨ ur ω? t(i, j, kt) vom Rang kt,max +1. In der diskretisierten Wassergehaltsgleichung (6.62) treten Terme auf, die einen hochgestellten, eingeklammerten Stern tragen (z. B. H(?) t). Hierbei handelt es sich um provisorische Werte f¨ ur den neuen Zeitpunkt des zeitlichen Iterationsschritts, die im Rahmen der probeweisen CTS-Verschiebungen zur Bestimmung von deren Position (siehe n¨ achster Abschnitt) berechnet werden; die endg¨ ultigen neuen Werte werden erst bei der Topographie-Berechnung bestimmt. Der Wasserdrainage-Funktion D(ω) ist in (6.62) noch nicht Rechnung getragen. Sie wird gem¨ aß ihrer Darstellung 99 Gleichungssysteme bei der L¨ osung auf. Mit impliziten x-Ableitungen erh¨ alt man h?(i, j)−h(i, j) ∆t =1 ∆x ¯ D? H(i+1 2, j)h?(i+ 1, j)−h?(i, j) ∆x−¯ D? H(i−1 2, j)h?(i, j)−h?(i−1, j) ∆x! +1 ∆y ¯ D? H(i, j +1 2)h(i, j + 1) −h(i, j) ∆y−¯ D? H(i, j −1 2)h(i, j)−h(i, j −1) ∆y! +(a⊥ s)?(i, j) + ∂b ∂t ? (i, j)−Pw b(i, j) ρ,(6.91) was f¨ ur jedes j= 1 . . . jmax −1 getrennt ein tridiagonales Gleichungssystem in den Unbekannten h?(i, j), (i= 0 . . . imax) darstellt (analog f¨ ur implizite y-Ableitungen, wird nicht extra gezeigt). An den Berandungen i= 0, imax bzw. j= 0, jmax des Rechengebietes wird bei der L¨ osung h?(i, j) = b?(i, j) vorausgesetzt, d. h., das Rechengebiet muß das simulierte Eisschild ganz umfassen. Mit den gewonnenen b?(i, j) und h?(i, j) und ggfs. den aus der Positionierungsroutine f¨ ur die CTS (n¨ achster Abschnitt) erhaltenen provisorischen H(?) t(i, j) k¨ onnen nun die neuen Dicken H?(i, j), H? c(i, j) und H? t(i, j) sowie die neue CTS-Position z? m(i, j) berechnet werden: Fall CB, TB: H?(i, j) = H? c(i, j) = h?(i, j)−b?(i, j), H? t(i, j) = 0, z? m(i, j) = b?(i, j).(6.92) Fall TL: H?(i, j) = h?(i, j)−b?(i, j), H? t(i, j) = H(?) t(i, j)H?(i, j) H(i, j), z? m(i, j) = b?(i, j) + H? t(i, j), H? c(i, j) = h?(i, j)−z? m(i, j).(6.93) Hierbei wird die provisorische Dicke der temperierten Bereiche H(?) t(i, j), welche noch unter dem Einfluß der alten Topographie bestimmt wurde, im Verh¨ altnis der neuen zur alten Gesamtdicke gestreckt. Nachdem nun die gesamte neue Topographie bekannt ist, k¨ onnen die noch fehlenden Zeitableitungen von h,zm,Hund Hc(f¨ ur die Zeitableitungen von bund Htsiehe (6.86) bzw. (6.107)) sowie alle r¨ aumlichen Ableitungen von h,zm,b,H,Hcund Ht 106 berechnet werden: ∂h ∂t ? (i, j) = h?(i, j)−h(i, j) ∆t, ∂H ∂t ? (i, j) = ∂h ∂t ? (i, j)−∂b ∂t ? (i, j), ∂zm ∂t ? (i, j) = ∂b ∂t ? (i, j) + ∂Ht ∂t ? (i, j), ∂Hc ∂t ? (i, j) = ∂h ∂t ? (i, j)−∂zm ∂t ? (i, j),(6.94) sowie ∂h ∂x ? (i+1 2, j) = h?(i+ 1, j)−h?(i, j) ∆x, ∂h ∂y ? (i, j +1 2) = h?(i, j + 1) −h?(i, j) ∆y(6.95) (entsprechend f¨ ur die r¨ aumlichen Ableitungen von zm,b,H,Hc,Ht). Schließlich wird f¨ ur den Fall TL der Volumenfluß durch die CTS a⊥ mbestimmt. Aufgrund der Vorgehensweise bei der CTS-Positionierung wird dessen station¨ arer Anteil a⊥ m,st (ohne Einfluß der CTS-Eigenbewegung) gesondert berechnet, (a⊥ m,st)?(i, j) = ¯v? x,c(i, j, 0)∂zm ∂x ? (i, j) + ¯v? y,c(i, j, 0)∂zm ∂y ? (i, j)−¯v? z,c(i, j, 0),(6.96) das vollst¨ andige a⊥ mfolgt hiermit gem¨ aß (a⊥ m)?(i, j) = (a⊥ m,st)?(i, j) + ∂zm ∂t ? (i, j).(6.97) Schmelztemperatur Diagnostische Berechnung aus der soeben bestimmten Topographie: Kalter Eisbereich: T? M,c(i, j, kc) = T0−βH? c(i, j)1−E(ζc(kc)).(6.98) Temperierter Eisbereich: T? M,t(i, j, kt) = T0−βH? c(i, j) + H? t(i, j) (1 −ζt(kt)).(6.99) 107 Basale Schmelzrate Diagnostische Berechnung mit Hilfe der neuen Werte f¨ ur Temperatur, Geschwindigkeit und Topographie: Fall CB: (Pw b)?(i, j)=0.(6.100) Fall TB: (Pw b)?(i, j) = κ(T? c(i, j, 0)) (ea−1) LaH? c(i, j)·T? c(i, j, 1) −T? c(i, j, 0) ∆ζc −κr LHr·T? r(i, j, kr,max)−T? r(i, j, kr,max −1) ∆ζr −ρg L¯v? x,c(i, j, 0) H? c(i, j)∂h ∂x ? (i, j) −ρg L¯v? y,c(i, j, 0) H? c(i, j)∂h ∂y ? (i, j).(6.101) Fall TL: (Pw b)?(i, j) = βκ(T? M,t(i, j, 0)) L−κr LHr·T? r(i, j, kr,max)−T? r(i, j, kr,max −1) ∆ζr −ρg L¯v? x,t(i, j, 0) H?(i, j)∂h ∂x ? (i, j) −ρg L¯v? y,t(i, j, 0) H?(i, j)∂h ∂y ? (i, j); (6.102) hierzu wird die Gesamtmenge des aus der temperierten Schicht der Eiss¨ aule (i, j) drainierten Wassers (parameterisiert durch die Drainagefunktion D(ω)) addiert. 6.5 Positionierung der CTS Bei der Berechnung von Temperatur und ggfs. Wassergehalt f¨ ur eine bestimmte Eiss¨ aule (i, j) tritt das Problem auf zu bestimmen, ob zum neuen Zeitpunkt des Iterationsschritts Fall CB (kalte Eisbasis), Fall TB (temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht) oder Fall TL (temperierte Eisbasis mit ¨ uberlagerter temperierter Schicht) vorliegt, und wo sich im letzteren Fall die CTS befindet. Dies wird mit einem Try-and-Error-Verfahren durchgef¨ uhrt, wobei die Vorgehensweise unterschiedlich ist, je nachdem von welchem der drei F¨ alle zum alten Zeitpunkt der Iteration ausgegangen wird. 108 Ausgangspunkt Fall CB f¨ ur S¨ aule (i,j) Schritt 1: Testrechnung unter der Annahme des Falles CB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53) – (6.58). Weiter mit Schritt 2. Schritt 2: Abfrage, ob die erhaltene Eisbasistemperatur den Druckschmelzpunkt unterschreitet, T? c(i, j, 0) ? < T0−βH(i, j). Falls ja, Verfahren beendet. Falls nein, weiter mit Schritt 3. Schritt 3: Testrechnung unter der Annahme des Falles TB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53), (6.54), (6.57) – (6.60). Weiter mit Schritt 4. Schritt 4: Abfrage, ob die erhaltene vertikale Temperaturableitung an der Eisbasis kleiner ist als der Clausius-Clapeyron-Gradient, ea−1 aHc(i, j)·T? c(i, j, 1) −T? c(i, j, 0) ∆ζc ? < β. Falls ja, Verfahren beendet. Falls nein, weiter mit Schritt 5. Schritt 5: Testrechnung unter der Annahme des Falles TL zum neuen Zeitpunkt, wobei die CTS um Hoffset t:= 1 mm ¨ uber die Eisbasis gelegt wird. Es werden also die provisorischen Topographiegr¨ oßen zu H(?) t(i, j) = Hoffset t, z(?) m(i, j) = b(i, j) + Hoffset t, H(?) c(i, j) = H(i, j)−Hoffset t gesetzt und das Problem (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66) gel¨ ost. Weiter mit Schritt 6. Schritt 6: Abfrage, ob die erhaltene Temperatur an der Unterseite des Kalteisbereichs den Druckschmelzpunkt unterschreitet, T? c(i, j, 0) ? < T0−βH(?) c(i, j). Falls ja, Verfahren beendet. Falls nein, weiter mit Schritt 7. 109 Schritt 7: Iterative Verschiebung der CTS um jeweils zshift m:= 1 m nach oben gem¨ aß H(?) t(i, j) = H(?) t(i, j) + zshift m, z(?) m(i, j) = z(?) m(i, j) + zshift m, H(?) c(i, j) = H(?) c(i, j)−zshift m, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), hiermit L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66); diese Verschiebungs-Prozedur wird gestoppt, wenn entweder die CTS um weniger als zshift mvon der Eisoberfl¨ ache entfernt ist, z(?) m(i, j)? > h(i, j)−zshift m, oder die erhaltene Temperatur an der Unterseite des Kalteisbereichs den Druckschmelzpunkt unterschreitet, T? c(i, j, 0) ? < T0−βH(?) c(i, j). Bei Erf¨ ullung des ersten Kriteriums wird das Verfahren beendet, bei Erf¨ ullung des zweiten weiter mit Schritt 8. Schritt 8: Interpolation der endg¨ ultigen CTS-Position aus der letzten (za:= z(?) m(i, j)) und vorletzten (zb:= z(?) m(i, j)−zshift m) provisorischen CTS-Position der Verschiebe-Prozedur, gewichtet mit den zugeh¨ origen Differenzen der errechneten Temperaturen an der Unterseite des Kalteisbereichs vom Druckschmelzpunkt, diffa,diffb := T? c(i, j)−(T0−βH(?) c(i, j)) (f¨ ur die CTS-Positionen zabzw. zb). Hierbei gilt diffa <0 und diffb >0; die 110 Interpolationsvorschrift lautet interpol := diffa diffb −diffa ·zshift m(<0), H(?) t(i, j) = H(?) t(i, j) + interpol, z(?) m(i, j) = z(?) m(i, j) + interpol, H(?) c(i, j) = H(?) c(i, j)−interpol, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), hiermit L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66). Beendigung des Verfahrens. Ausgangspunkt Fall TB f¨ ur S¨ aule (i,j) Schritt 1: Testrechnung unter der Annahme des Falles TB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53), (6.54), (6.57) – (6.60). Weiter mit Schritt 2. Schritt 2: Abfrage, ob die erhaltene vertikale Temperaturableitung an der Eisbasis kleiner ist als der Clausius-Clapeyron-Gradient, ea−1 aHc(i, j)·T? c(i, j, 1) −T? c(i, j, 0) ∆ζc ? < β. Falls ja, weiter mit Schritt 3a. Falls nein, weiter mit Schritt 3b. Schritt 3a: Testrechnung unter der Annahme des Falles CB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53) – (6.58). Weiter mit Schritt 4a. Schritt 4a: Abfrage, ob die erhaltene Basistemperatur den Druckschmelzpunkt unterschreitet, T? c(i, j, 0) ? < T0−βH(i, j). Falls ja, Verfahren beendet. Falls nein, weiter mit Schritt 5a. 111 Schritt 5a: Rechnung unter der Annahme des Falles TB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53), (6.54), (6.57) – (6.60); Beendigung des Verfahrens. Schritt 3b: Annahme des Falles TL zum neuen Zeitpunkt, Durchf¨ uhrung der Schritte 5-8 von “Ausgangspunkt Fall CB f¨ ur S¨ aule (i, j)” zur Positionierung der CTS, Beendigung des Verfahrens. Ausgangspunkt Fall TL f¨ ur S¨ aule (i,j) Schritt 1: Testrechnung unter der Annahme des Falles TL zum neuen Zeitpunkt mit der alten CTS-Position, d. h., L¨ osung von (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66). Weiter mit Schritt 2. Schritt 2: Abfrage, ob die erhaltene Temperatur an der Unterseite des Kalteisbereichs den Druckschmelzpunkt ¨ uberschreitet, T? c(i, j, 0) ? > T0−βHc(i, j). Falls ja, weiter mit Schritt 3a. Falls nein, weiter mit Schritt 3b. Schritt 3a: Durchf¨ uhrung der Schritte 5-8 von “Ausgangspunkt Fall CB f¨ ur S¨ aule (i, j)” zur Positionierung der CTS (Verschiebung nach oben), Beendigung des Verfahrens. Schritt 3b: Iterative Verschiebung der CTS um jeweils zshift m:= 1 m nach unten gem¨ aß H(?) t(i, j) = H(?) t(i, j)−zshift m, z(?) m(i, j) = z(?) m(i, j)−zshift m, H(?) c(i, j) = H(?) c(i, j) + zshift m, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), hiermit L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66); 112 diese Verschiebungs-Prozedur wird gestoppt, wenn entweder die CTS um weniger als zshift mvon der Eisbasis entfernt ist, z(?) m(i, j)? < b(i, j) + zshift m, oder die erhaltene Temperatur an der Unterseite des Kalteisbereichs den Druckschmelzpunkt ¨ uberschreitet, T? c(i, j, 0) ? > T0−βH(?) c(i, j). Bei Erf¨ ullung des ersten Kriteriums weiter mit Schritt 4bα, bei Erf¨ ullung des zweiten weiter mit Schritt 4bβ. Schritt 4bα:Positionierung der CTS um Hoffset t:= 1 mm ¨ uber der Eisbasis, H(?) t(i, j) = Hoffset t, z(?) m(i, j) = b(i, j) + Hoffset t, H(?) c(i, j) = H(i, j)−Hoffset t, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66). Weiter mit Schritt 5bα. Schritt 5bα:Abfrage, ob die erhaltene Temperatur an der Unterseite des Kalteisbereichs den Druckschmelzpunkt unterschreitet, T? c(i, j, 0) ? < T0−βH(?) c(i, j). Falls ja, weiter mit Schritt 6bα.1. Falls nein, weiter mit Schritt 6bα.2. Schritt 6bα.1: Rechnung unter der Annahme des Falles TB zum neuen Zeitpunkt, d. h., L¨ osung von (6.53), (6.54), (6.57) – (6.60). Beendigung des Verfahrens. 113 Schritt 6bα.2: Interpolation der endg¨ ultigen CTS-Position aus der letzten (za:= z(?) m(i, j) = b(i, j) + Hoffset t) und vorletzten (zb) provisorischen CTS-Position der Verschiebe-Prozedur, gewichtet mit den zugeh¨ origen Differenzen der errechneten Temperaturen an der Unterseite des Kalteisbereichs vom Druckschmelzpunkt, diffa,diffb := T? c(i, j)−(T0−βH(?) c(i, j)) (f¨ ur die CTS-Positionen zabzw. zb). Hierbei gilt diffa >0 und diffb <0; die Interpolationsvorschrift lautet interpol := diffa diffa −diffb ·(zb−za) (>0), H(?) t(i, j) = H(?) t(i, j) + interpol, z(?) m(i, j) = z(?) m(i, j) + interpol, H(?) c(i, j) = H(?) c(i, j)−interpol, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), hiermit L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66). Beendigung des Verfahrens. Schritt 4bβ:Interpolation der endg¨ ultigen CTS-Position aus der letzten (za:= z(?) m(i, j)) und vorletzten (zb:= z(?) m(i, j)+zshift m) provisorischen CTS-Position der Verschiebe-Prozedur, gewichtet mit den zugeh¨ origen Differenzen der errechneten Temperaturen an der Unterseite des Kalteisbereichs vom Druckschmelzpunkt, diffa,diffb := T? c(i, j)−(T0−βH(?) c(i, j)) (f¨ ur die CTS-Positionen zabzw. zb). Hierbei gilt diffa >0 und diffb <0; die 114 Interpolationsvorschrift lautet interpol := diffa diffa −diffb ·zshift m(>0), H(?) t(i, j) = H(?) t(i, j) + interpol, z(?) m(i, j) = z(?) m(i, j) + interpol, H(?) c(i, j) = H(?) c(i, j)−interpol, ∂Ht ∂t (?) (i, j) = H(?) t(i, j)−Ht(i, j) f ∆t, ∂zm ∂t (?) (i, j) = ∂b ∂t(i, j) + ∂Ht ∂t (?) (i, j), ∂Hc ∂t (?) (i, j) = ∂h ∂t (i, j)−∂zm ∂t (?) (i, j), (a⊥ m)(?)(i, j) = a⊥ m,st(i, j) + ∂zm ∂t (?) (i, j), hiermit L¨ osung des Fall-TL-Problems (6.53), (6.54), (6.57), (6.58), (6.60) – (6.66). Beendigung des Verfahrens. Gl¨ attung der CTS Die Durchf¨ uhrung des oben beschriebenen Verfahrens zur Unterscheidung zwischen Fall CB, Fall TB und Fall TL sowie bei letzterem zur Positionierung der CTS f¨ ur alle Eiss¨ aulen (i, j) des Rechengebietes tendiert dazu, r¨ aumliche Oszillationen des CTSVerlaufes von Gitterpunkt zu Gitterpunkt zu produzieren. Daher wird nach Abarbeitung aller Eiss¨ aulen eine numerische Gl¨ attung der erhaltenen Dicken der temperierten Bereiche H(?) t(i, j) vorgenommen, wobei auf die vier n¨ achsten Nachbarn zugegriffen wird. Im ersten Schritt wird hierzu das Volumen temperierten Eises Vtemp vor der Gl¨ attung berechnet, Vtemp = imax−1 X i=1 jmax−1 X j=1 H(?) t(i, j) ∆x∆y. (6.103) Der zweite Schritt besteht in der eigentlichen Gl¨ attung, welche jedoch nur unter der Voraussetzung einer temperierten Basis zum neuen Zeitpunkt ausgef¨ uhrt wird, H(?) t,smooth(i, j) = (1 −4DHt)H(?) t(i, j) +DHtH(?) t(i+ 1, j) + DHtH(?) t(i−1, j) +DHtH(?) t(i, j + 1) + DHtH(?) t(i, j −1) (Fall TB, TL), H(?) t,smooth(i, j) = 0 (Fall CB).(6.104) 115 •Oberfl¨ achentemperatur-Antrieb ∆Ts, basale homologe Temperatur T0 b,div bei (xdiv, ydiv), gesamtes Eisvolumen Vges, maximale Eish¨ ohe hmax, temperiertes Eisvolumen Vtemp, maximale Dicke der temperierten Eisschicht Ht,max, eisbedeckte Bodenfl¨ ache Ai,b, von temperiertem Eis bedeckte Bodenfl¨ ache At,b, Massenfluß qx,half bei (xhalf , ydiv), jeweils als Funktion der Zeit t. Hierbei bedeuten xdiv := 750 km und ydiv := 750 km die Position des Scheitelpunktes des Eisschildes (Mittelpunkt des Quadrates) sowie xhalf := 1150 km eine Position auf halbem Weg von xdiv bis zum rechten Eisrand. Abgesehen von den ZeitentwicklungsPlots beziehen sich alle Gr¨ oßen auf den erreichten Endzustand der betreffenden Simulation. Im Vergleich der beiden Steady-State-L¨ aufe ssfml1 und ssfml2 zeigen sich einige charakteristische Unterschiede. So erkennt man beim Vergleich der resultierenden Topographien (Abb. 7.1, 7.4, 7.7, 7.10), daß zwar beide L¨ aufe ein kappenf¨ ormiges Eisschild produzieren, die Eiskappe des Level 2-Laufes jedoch wesentlich spitzer erscheint. Das liegt daran, daß bei Level 2 basales Gleiten auf temperierter Eisbasis ber¨ ucksichtigt ist, wodurch das randnahe Eis schneller zu den Seiten hin abfließt als bei Level 1 und die Randbereiche somit d¨ unner werden. Damit einher gehen auch die bei ssfml2 gr¨ oßeren Oberfl¨ achengeschwindigkeiten (Abb. 7.2, 7.5, 7.8, 7.11). Eine zus¨ atzliche Konsequenz des basalen Gleitens bei Level 2 ist die durch die generell gr¨ oßeren Geschwindigkeiten im Eis verst¨ arkte Advektion von kaltem Oberfl¨ acheneis nach unten, wodurch ssfml2 ein k¨ alteres Eisschild als ssfml1 hervorbringt (Abb. 7.3, 7.4, 7.5, 7.6, 7.9, 7.10, 7.11, 7.12). Sehr kraß tritt dieser Unterschied bei der Ausdehnung der temperierten Bereiche zutage. Zwar unterscheiden sich die Bodenfl¨ achen mit temperierter Eisbasis At,b mehr in der Geometrie als in der Gr¨ oße (Abb. 7.3, 7.9), die temperierten Eisvolumina sind jedoch etwa um den Faktor 20 verschieden (Abb. 7.4, 7.6, 7.10, 7.12). Dies zeigt, daß eine gute Kenntnis der physikalischen Vorg¨ ange an der Eisbasis sehr wichtig f¨ ur die Vorhersage der Ausdehnung einer temperierten Eisschicht ist. Ein weiterer bemerkenswerter Unterschied zwischen ssfml1 und ssfml2 besteht im Verlauf des Geschwindigkeitsprofils vx,half (Abb. 7.5, 7.11). Hier ist der Grenzschichtcharakter des bodennahen, großen Scherkr¨ aften unterliegenden Eises f¨ ur den Level 2Lauf wesentlich st¨ arker ausgepr¨ agt als beim Level 1-Lauf, bei dem der Anstieg der Geschwindigkeit mit der H¨ ohe deutlich gleichm¨ aßiger verl¨ auft. Das hat seine Ursache darin, daß bei Level 2 der starke Anstieg der Eisfluidit¨ at mit der Temperatur ber¨ ucksichtigt ist, wodurch die Scherung des relativ warmen, bodennahen Eises gegen¨ uber dem k¨ alteren Eis weiter oben beg¨ unstigt wird. Der Verlauf des Wassergehaltes ω2scheint bei ssfml1 (Abb. 7.5) eine Inkonsi122 stenz aufzuweisen. Obwohl bei der zugeh¨ origen Position (x2, ydiv) Schmelzbedingungen an der CTS herrschen (Abb. 7.4), geht der Wassergehalt an der CTS nicht auf Null zur¨ uck. Die Ursache f¨ ur dieses Verhalten liegt in der Gl¨ attungsprozedur f¨ ur die CTS-Position (§6.5) zusammen mit den bei Level 1 auftretenden sehr großen Dicken der temperierten Eisschicht aufgrund des unber¨ ucksichtigten basalen Gleitens (siehe oben). Die Gl¨ attung nach erfolgter CTS-Positionierung bewirkt, daß bei (x2, ydiv) die relativ kleine Dicke Htder temperierten Eisschicht durch gr¨ oßere Dicken in der Nachbarschaft vergr¨ oßert wird. In der darauffolgenden Iteration hat die CTS daher die Tendenz, wieder nach unten zu wandern, wodurch bei der Berechnung des Wassergehaltes Gefrierbedingungen diagnostiziert werden. Dieses Ph¨ anomen tritt jedoch bei der physikalisch realistischeren Level 2-Simulation nicht auf, weil dort wegen der wesentlich geringeren Dicken der temperierten Eisschicht der absolute Einfluß der CTS-Gl¨ attung auf Htkleiner ist. In Abb. 7.11 erkennt man daher auch einen an der CTS auf Null zur¨ uckgehenden Wassergehalt ω2, wie es den dort vorherrschenden Schmelzbedingungen entspricht. Bei den vier transienten L¨ aufen t2fml1, t4fml1, t2fml2, t4fml2 erkennt man zun¨ achst, daß in der zeitlichen Entwicklung mit Ausnahme der sowieso konstanten Vereisungsfl¨ ache Ai,b alle gezeigten Gr¨ oßen der antreibenden Milankovi´c-Periode tMil = 20 bzw. 40 ka folgen (Abb. 7.15, 7.18, 7.21, 7.24). Das antreibende Sinussignal erscheint zwar bei einigen Gr¨ oßen verzerrt, was auf Oberschwingungen hindeutet, jedoch treten keine Effekte wie Periodenverdopplungen oder gar chaotisches Verhalten auf, wie dies von Abe-Ouchi [1] bei der Modellierung der EGIG-Traverse Gr¨ onlands gefunden wurde5. Es f¨ allt auf, daß sowohl f¨ ur Level 1 als auch f¨ ur Level 2 die Schwingungsamplituden der rein dynamischen Gr¨ oßen Vges,hmax und qx,half im Vergleich der beiden verschiedenen Anregungsperioden fast gleich sind, w¨ ahrend die Amplituden der direkt mit dem Temperaturfeld zusammenh¨ angenden Gr¨ oßen T0 b,div,Vtemp,Ht,max und At,b bei 40 ka-Anregung deutlich gr¨ oßer sind als bei 20 ka-Anregung. Das liegt daran, daß sich die Eistopographie wesentlich schneller auf ge¨ anderte Randbedingungen einstellt als das Temperaturfeld (vgl. auch die Zeitentwicklung bei den Steady-State-L¨ aufen; Abb. 7.6, 7.12), wodurch das Temperaturfeld der kleineren Periode deutlich schlechter folgen kann als der gr¨ oßeren. Es mag erstaunen, daß in allen vier F¨ allen das Eisvolumen Vges und die Maximalh¨ ohe hmax in Phase statt in Gegenphase mit der antreibenden Temperatur ∆Ts sind, w¨ urde man doch bei w¨ armeren Lufttemperaturen ein d¨ unneres Eisschild erwarten. Der Grund daf¨ ur ist in der speziellen Vorgabe der Randbedingungen zu su5Das bedeutet jedoch nicht, daß die einfache Geometrie des EISMINT-Eisschildes prinzipiell kein irregul¨ ares Verhalten zul¨ aßt. Es besagt lediglich, daß dies im Rahmen der hier gew¨ ahlten transienten klimatischen Randbedingungen nicht auftritt. 123 chen, denn erstens ist kein Oberfl¨ achenschmelzen vorgegeben, welches bei w¨ armeren Temperaturen verst¨ arkt Eis abbauen w¨ urde, und zweitens nimmt gem¨ aß (7.2) die Akkumulationsrate mit steigender Oberfl¨ achentemperatur stark zu. Diese Zunahme ¨ uberkompensiert offenbar die leichtere Eisdeformierbarkeit aufgrund der h¨ oheren Temperaturen und des gr¨ oßeren temperierten Eisvolumens. Vergleicht man die Endzust¨ ande der transienten L¨ aufe mit denen der SteadyState-L¨ aufe, so zeigt sich, daß bei den transienten L¨ aufen ausgepr¨ agte oberfl¨ achennahe Temperaturinversionen (d. h., die Temperatur nimmt mit steigender Tiefe ab statt zu) sowohl im Randbereich als auch in der Mitte des Eisschildes auftreten (Abb. 7.13, 7.16, 7.19, 7.22), w¨ ahrend dies bei den Steady-State-L¨ aufen nur f¨ ur den Randbereich gilt (Abb. 7.4, 7.10). Besonders augenf¨ allig wird das durch die geschlossenen −35◦CIsolinien der homologen Temperatur. Diese Erscheinung liegt darin begr¨ undet, daß die letzte Viertelphase der Anregung vor dem Ende des jeweiligen Laufes nach t= 200 ka in allen F¨ allen einen Temperaturanstieg von ∆Ts=−10◦C auf ∆Ts= 0◦C beinhaltet (welcher im ¨ ubrigen den ¨ Ubergang von der letzten (W¨ urm-) Kaltzeit zur heutigen (Neo-) Warmzeit ziemlich realistisch wiedergibt). Das Eis in einiger Tiefe hat diese niedrigen Temperaturen noch gespeichert, woraus die beobachtete Inversion resultiert. Bei den Steady-State-L¨ aufen l¨ aßt sich die nahe dem Eisrand ebenfalls vorhandene Temperaturinversion nat¨ urlich nicht so erkl¨ aren; hier liegt die Ursache in der Form der Eisbewegung selbst, welche kaltes Oberfl¨ acheneis aus den randferneren Gebieten advektiv unter das w¨ armere Oberfl¨ acheneis in Randn¨ ahe bef¨ ordert (selbstverst¨ andlich tritt dieser Effekt auch bei den transienten L¨ aufen auf und tr¨ agt dort ebenfalls zur Temperaturinversion in Randn¨ ahe bei). Die bereits oben bei der Besprechung der Steady-State-L¨ aufe erw¨ ahnte Erscheinung, daß der Wassergehalt ω2bei ssfml1 trotz Schmelzbedingungen an der CTS nicht verschwindet, tritt auch bei t2fml1 und t4fml1 auf (Abb. 7.14, 7.17) und hat die gleiche Ursache. Bei den realistischeren Level 2-L¨ aufen t2fml2 und t4fml2 (Abb. 7.20, 7.23) geht hingegen ω2an der CTS auf Null zur¨ uck; t4fml2 ergibt zudem als einziger der sechs betrachteten L¨ aufe ein Wassergehaltsprofil ω2, welches auf einigen Gitterpunkten auch Werte zwischen Null und ωmax = 1% annimmt. Abschließend sei auf eine deutlich sichtbare Auswirkung der bei den transienten Level 2-L¨ aufen ber¨ ucksichtigten thermischen Tr¨ agheit der Lithosph¨ are hingewiesen. Betrachtet man die Zeitentwicklung der homologen Basaltemperatur unter dem Scheitelpunkt des Eisschildes T0 b,div, so sieht man, daß diese bei den beiden Level 1-L¨ aufen (Abb. 7.15, 7.18) bereits nach ca. 75 ka gut eingeschwungen ist, w¨ ahrend f¨ ur die Level 2-L¨ aufe (Abb. 7.21, 7.24) selbst am Ende der jeweiligen Simulation nach 200 ka noch eine Drift zu gr¨ oßeren Werten hin festzustellen ist. Solche extrem großen Zeit124 konstanten bei der Einstellung der Temperatur auf ge¨ anderte Randbedingungen sind typisch f¨ ur die Auswirkung der Lithosph¨ are und zeigen auf, daß diese f¨ ur realistische Simulationen transienter Temperaturentwicklungen in Eisschilden (und damit einhergehend auch temperierter Eisbereiche) unbedingt ber¨ ucksichtigt werden muß. 125 126 7.3 Abbildungen zu den Simulationen f¨ ur das EISMINTEisschild 127 128 Abbildung 7.1: Endzustand von Lauf ssfml1: Topographie der Eisoberfl¨ ache (in km ¨ uber Meeresh¨ ohe). Die H¨ ohendifferenz zwischen den Isolinien betr¨ agt 200 m. 129 Abbildung 7.2: Endzustand von Lauf ssfml1: Eisoberfl¨ achengeschwindigkeit (in km/a). Die Geschwindigkeiten zu aufeinanderfolgenden Isolinien unterscheiden sich jeweils um den Faktor zwei. 130 Abbildung 7.3: Endzustand von Lauf ssfml1: Homologe Temperatur an der Eisbasis (in ◦C). Die Temperaturdifferenz zwischen den Isolinien betr¨ agt 3◦C. Bei den mit Symbolen markierten Punkten befindet sich die Eisbasis auf dem Druckschmelzpunkt, wobei leere Diamantsymbole eine temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht, ausgef¨ ullte Diamantsymbole eine temperierte Schicht mit Schmelzbedingungen an der CTS, ausgef¨ ullte Kreise eine temperierte Schicht mit Gefrierbedingungen an der CTS bedeuten. 131 Abbildung 7.10: Endzustand von Lauf ssfml2: Querschnitt bei y= 750 km. Oben: Eisgeschwindigkeit. Mitte: homologe Eistemperatur (in ◦C). Unten: Dicke der basalen temperierten Eisschicht (leere Kreise: kalte Eisbasis; leere Diamantsymbole: temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht; ausgef¨ ullte Diamantsymbole: temperierte Schicht mit Schmelzbedingungen an der CTS; ausgef¨ ullte Kreise: temperierte Schicht mit Gefrierbedingungen an der CTS.) 138 Abbildung 7.11: Endzustand von Lauf ssfml2: Querschnitt bei y= 750 km. vx,s,T0 b,qxals Funktion von x;vx,half ,ω2(mit Angabe der zugeh¨ origen Position x2), vz,div,vz,half (gestrichelt), T0 div,T0 half als Funktion von z(f¨ ur ω2auf Ht, f¨ ur die anderen Plots auf Hnormiert). Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 139 Abbildung 7.12: Lauf ssfml2: Zeitliche Entwicklung von ∆Ts,T0 b,div,Vges,hmax,Vtemp,Ht,max, Ai,b,At,b (gestrichelt), qx,half . Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 140 Abbildung 7.13: Endzustand von Lauf t2fml1: Querschnitt bei y= 750 km. Oben: Eisgeschwindigkeit. Mitte: homologe Eistemperatur (in ◦C). Unten: Dicke der basalen temperierten Eisschicht (leere Kreise: kalte Eisbasis; leere Diamantsymbole: temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht; ausgef¨ ullte Diamantsymbole: temperierte Schicht mit Schmelzbedingungen an der CTS; ausgef¨ ullte Kreise: temperierte Schicht mit Gefrierbedingungen an der CTS.) 141 Abbildung 7.14: Endzustand von Lauf t2fml1: Querschnitt bei y= 750 km. vx,s,T0 b,qxals Funktion von x;vx,half ,ω2(mit Angabe der zugeh¨ origen Position x2), vz,div,vz,half (gestrichelt), T0 div, T0 half als Funktion von z(f¨ ur ω2auf Ht, f¨ ur die anderen Plots auf Hnormiert). Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 142 Abbildung 7.15: Lauf t2fml1: Zeitliche Entwicklung von ∆Ts,T0 b,div,Vges,hmax,Vtemp,Ht,max, Ai,b,At,b (gestrichelt), qx,half . Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 143 Abbildung 7.16: Endzustand von Lauf t4fml1: Querschnitt bei y= 750 km. Oben: Eisgeschwindigkeit. Mitte: homologe Eistemperatur (in ◦C). Unten: Dicke der basalen temperierten Eisschicht (leere Kreise: kalte Eisbasis; leere Diamantsymbole: temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht; ausgef¨ ullte Diamantsymbole: temperierte Schicht mit Schmelzbedingungen an der CTS; ausgef¨ ullte Kreise: temperierte Schicht mit Gefrierbedingungen an der CTS.) 144 Abbildung 7.17: Endzustand von Lauf t4fml1: Querschnitt bei y= 750 km. vx,s,T0 b,qxals Funktion von x;vx,half ,ω2(mit Angabe der zugeh¨ origen Position x2), vz,div,vz,half (gestrichelt), T0 div, T0 half als Funktion von z(f¨ ur ω2auf Ht, f¨ ur die anderen Plots auf Hnormiert). Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 145 Abbildung 7.18: Lauf t4fml1: Zeitliche Entwicklung von ∆Ts,T0 b,div,Vges,hmax,Vtemp,Ht,max, Ai,b,At,b (gestrichelt), qx,half . Die Bedeutung der einzelnen Gr¨ oßen ist im Haupttext erkl¨ art. 146 Abbildung 7.19: Endzustand von Lauf t2fml2: Querschnitt bei y= 750 km. Oben: Eisgeschwindigkeit. Mitte: homologe Eistemperatur (in ◦C). Unten: Dicke der basalen temperierten Eisschicht (leere Kreise: kalte Eisbasis; leere Diamantsymbole: temperierte Eisbasis ohne ¨ uberlagerte temperierte Schicht; ausgef¨ ullte Diamantsymbole: temperierte Schicht mit Schmelzbedingungen an der CTS; ausgef¨ ullte Kreise: temperierte Schicht mit Gefrierbedingungen an der CTS.) 147