Die Verwendung des footprint-Modells nach SCHMID (1997) zur stabilitätsabhängigen Bestimmung der Rauhigkeitslänge

UNIVERSITÄT BAYREUTH Abt. Mikrometeorologie Die Verwendung des footprint-Modells nach SCHMID (1997) zur stabilitätsabhängigen Bestimmung der Rauhigkeitslänge Mathias Göckede Arbeitsergebnisse Nr.15 Bayreuth, März 2001 Arbeitsergebnisse, Universität Bayreuth, Abt. Mikrometeorologie, Print, ISSN 1614-8916 Arbeitsergebnisse, Universität Bayreuth, Abt. Mikrometeorologie, Internet, ISSN 1614-8924 http://www.bayceer.uni-bayreuth.de/mm/ Eigenverlag: Universität Bayreuth, Abt. Mikrometeorologie Vervielfältigung: Druckerei der Universität Bayreuth Herausgeber: Prof. Dr. Thomas Foken Universität Bayreuth, Abteilung Mikrometeorologie 0-95440 Bayreuth Die Verantwortung über den Inhalt liegt beim Autor. Die Verwendung des footprint-Modells nach SCHMID (1997) zur stabilitätsabhängigen Bestimmung der Rauhigkeitslänge 1 Einleitung Bei der Modellierung von Austauschströmen in der atmosphärischen Bodenschicht wird für die Berechnung meist vorausgesetzt, dass die Oberfläche weitgehend homogen und mehr oder weniger flach ist. Die meisten meteorologischen Experimente werden daher an Standorten durchgeführt, die diese Eigenschaften erfüllen, denn in diesem Fall können die Charakteristika der Turbulenz gut nachvollzogen werden. Will man die unter solchen idealen Bedingungen gewonnenen Parametrisierungsansätze übertragen auf komplexer gegliederte Umgebungen, muss auch deren räumliche Heterogenität berücksichtigt werden. Im typischen ländlichen Raum liegen beispielsweise Flächen unterschiedlicher Landnutzung in einer Art Flickenteppichmuster nebeneinander. Dabei haben Wälder völlig andere Einflüsse auf das Windfeld als Ackerflächen oder Weiden, was sich in voneinander abweichenden Rauhigkeitslängen wiederspiegelt. Zudem kann sich jede Landnutzungsart anders auf den fühlbaren Wärmestrom auswirken, da in Abhängigkeit sowohl von der Vegetationsart und -zusammensetzung als auch von den Feuchtebedingungen des Untergrundes unterschiedliche Quellen bzw. Senken für Wärme entstehen. Generell tendieren die Turbulenz-Charakteristika in der atmosphärischen Bodenschicht immer gegen ein Gleichgewicht mit den Bedingungen der Unterlage. Ist das Gelände aber nicht einheitlich, sondern weist Veränderungen in seiner generellen Gestalt auf, können die Messungen nicht auf die Eigenschaften des Standorts selbst bezogen werden. Denn unter dem Einfluss von horizontaler Advektion und vertikaler Diffusion werden meteorologische Messungen nicht nur von einem Punkt an der Oberfläche beeinflusst, sondern von einer möglicherweise sehr großen Region, dem Quellgebiet. Demnach repräsentiert eine solche Messung das räumliche Mittel vieler Einzelquellen eines weiter entfernten Bereichs entgegen der Windrichtung. Die Abgrenzung und Gewichtung dieses Einflussgebiets der Unterlage auf einen Sensor, welcher sich in einer bestimmten Höhe über Grund an der betreffenden Position befindet, ist das Untersuchungsziel der sogenannten ,footprint-Analysen'. Das berechnete Areal wird als Quellgebiet der Bodeneigenschaften für einen spezifischen Messpunkt, im Englischen ,source-area', bezeichnet. Bezüglich des fühlbaren Wärmestroms beispielsweise enthält diese sourcearea alle Quellen oder Senken für Wärme, welche die Luftpakete beeinflussen, die das Messgerät passieren. Lage und Dimension des Gebiets verändern sich ständig in Abhängigkeit von Windrichtung und -geschwindigkeit sowie den Stabilitätsbedingungen der atmosphärischen Schichtung. Sind diese Parameter bekannt, so kann man anhand einer geeigneten Gewichtungsfunktion den innerhalb liegenden Flächen ihren jeweiligen Einfluss am Gesamteffekt zuordnen. Die Berechnung des in Frage kommenden Parameters, beispielsweise des fühlbaren Wärmestroms, erfolgt durch Integration über diese Funktion. Sie wird nach SCHMID (1994) einheitlich als ,source weight function' oder Quelldichtefunktion bezeichnet, nachdem vor dieser festen Definition auch der synonyme Ausdruck ,footprint function' verwendet worden war. Mit Hilfe der footprint-Analysen können also die Regionen der Bodenoberfläche identifiziert und in ihrem Einfluss bewertet werden, die ein Sensor tatsächlich im aerodynamischen Sinne ,sieht'. Voraussetzung für diese Theorie ist nach SCHMID & OKE (1990) allerdings, dass das Strömungsfeld auf 2 seinem Weg von der source-area zum Messstandort nicht nachhaltig gestört worden ist. Sowohl der mittlere Fluss als auch die turbulente Durchmischung sollten daher stationär und homogen sein. Da diese Bedingung über Regionen mit Vegetationsoder Rauhigkeitswechseln kaum realisierbar ist, muss vor jeder Anwendung der footprint-Modelle überprüft werden, ob das betreffende Gebiet die geforderten Anforderungen erfüllt, und inwieweit sich Störungen auf die Verlässlichkeit der Resultate auswirken. Des weiteren können nur Parameter modelliert werden, die sich passiv in einem homogenen Diffusionsfeld ausbreiten. Die source-area bezieht sich nur auf ein Gebiet am Erdboden, welches die turbulenten Austauschprozesse beeinflusst, Produktion oder Vernichtung in der Atmosphäre selbst werden nicht mit berücksichtigt. In der vorliegenden Untersuchung wurde ein footprint-Modell dazu verwendet, die effektive Rauhigkeitslänge für einen Sensor in beliebiger Höhe unter verschiedenen meteorologischen Bedingungen zu berechnen. Da das Quellgebiet fiir die Messgeräte an unterschiedlichen Höhen eines Mastes in seinen Dimensionen variiert, können sich in heterogenen Umgebungen auch die ermittelten Rauhigkeitswerte stark voneinander unterscheiden. Die individuelle Bestimmung der Bodenrauhigkeits-Einflüsse für jede Höhe soll beispielsweise ermöglichen, ein Windprofil zu berechnen, für welches die komplexen Auswirkungen einer uneinheitlichen Unterlage berücksichtigt worden sind. Unterschiedliche Varianten von footprint-Analysen wurden in den vergangenen Jahren von mehreren Autoren vorgestellt. Zu erwähnen sind hier die Ansätze nach SCHUEPP et al. (1990), LECLERC & THURTELL (1990), HORST & WEIL (1992, 1994) sowie SCHMID (1994, 1997). Der hier verwendete Berechnungsansatz stammt von SCHMID (1997). Er soll im folgenden Abschnitt näher vorgestellt werden. 2 Das footprint-Modell nach SCHMID (1997) Das im folgenden Abschnitt vorgestellte footprint-Modell wurde von SCHMID & OKE (1990) entwickelt, weitere Nachbesserungen fiihrten zu modifizierten Versionen nach SCHMID (1994) und SCHMID (1997). Es basiert auf dem Grundgedanken der Gauss'schen Diffusion, die einen passiven Parameter in einem horizontal homogenen atmosphärischen Strömungsfeld verteilt. Unter dieser Annahme emittiert eine künstliche Punktquelle eine Fahne in die unterste Luftschicht, welche sich im Verlauf des Transports immer mehr ausbreitet und dabei verdünnt. Platziert man einen Sensor in die Strömungsrichtung dieser Fahne, so wird sich die gemessene Konzentration je nach der relativen Position des Messgerätes zur Quelle ändern. Befindet sich dieser Sensor genau im Zentrum des Abluftstroms, ergeben sich die höchsten Messwerte für den Parameter. Die Quelle markiert damit die sogenannte ,maximum source location' (Punkt des maximalen Bodeneinflusses). Bewegt man den Ort der Emission von diesem Maximum in die Windrichtung oder quer dazu, so wird die ermittelte Konzentration immer mehr abnehmen, bis letztendlich eine Position erreicht wird, ab der die Auswirkungen der Punktquelle auf die Messung faktisch verschwindet. Die Summe aller dieser Stellen bilden eine geschlossene Kurve, welche als Grenze den äußeren Rand der source-area markiert. In meteorologischen Studien arbeitet man nun nicht nur mit einer einzelnen Quelle, deren Einfluss auf den Sensor abgeschätzt werden soll. Denn jeder einzelne Punkt innerhalb des Quellgebietes wirkt sich auf die Messung aus, beispielsweise als Wärmequelle oder -senke bei der Bestimmung des turbulenten fiihlbaren Wärmestroms. Allen Flächen innerhalb dieses Gebiets muss daher über eine geeignete Quelldichtefunktion ein gewichteter Anteil am Gesamteffekt auf den Sensor zugewiesen werden. Mathematisch wird dieses Verfahren umgesetzt, indem die Ausbreitungsfahne, die von der maximum source loeation ausgeht, bestimmt wird. Diese Fahne bildet am Messstandort in der festgelegten Höhe 3 des Sensors ein bestimmtes Verteilungsmuster der Konzentration, welches sich mit einem geeigneten Diffusionsmodell berechnen lässt. SCHMID (1997) verwendet dazu die sogenannte Dirac-DeltaVerteilungsfunktion. Nach SCHMID & OKE (1990) kann man diese Verteilungskurve nun übertragen auf eine künstliche Emissionsquelle am Fuße des Messmastes, welche von einem virtuellen Wind entgegen der eigentlichen Windrichtung ausgebreitet wird. Diese Ausbreitungswolke markiert direkt über dem Punkt des maximalen Einflusses eine Fläche, welche genau die source-area identifiziert, wenn man sie auf den Erdboden projiziert. Dieses Konzept ist in Abbildung 1 dargestellt. LI ~ Wind /"'T"\ I I \ " I \ -' I \ , , , _f,-------+--- .. , -- I ........... I I , x , " " , " " , , x t Height ,r.L__ ~>---;-,"'- ___ -= __ -- Zs , , , , , , " " , , iJ =-ü ..... -- Virtual wind " ----~~--+------~(r---4-~--------------~~.----------~ Ma)(. {oufee loeation ~--....v P - source area Mit: P = Anteil am Effekt des emittierten Parameters X = gemessener Effekt Xp = minimale messbare Effekt-Level Xma, = maximaler gemessener Effekt Zs = Höhe des Sensors über Grund [m 1 ··X-O--x Distance X' = gemessener Effekt der virtuellen Quelle OJ = Quelldichtefunktion für einen Einzelpunkt (Op = maximale source-area-Fläche n(x,y) = Quelldichte-Verteilungsfunktion Abbildung 1: Berechnung der source-area für einen Sensor durch eine reverse Ausbreitungsfahne gegen die Windrichtung. Darstellung als schematischer Querschnitt. Nach SCHMID & OKE (1990). Die Berechnung des Wertes eines zu bestimmenden Parameters erfolgt als Integration über alle Punktquellen im Anströmungsbereich des Messstandortes, modifiziert durch eine Wahrscheinlichkeits-Transferfunktion für die Diffusion. Nach SCHMID (1997) ergibt sich folgende generelle Beziehung: mit: 11 = Wert der zu bestimmenden Größe r = Position des Messmastes r' = Position der Punktquelle Q'l = Stärke der Punktquelle f(r-r') = Wahrscheinlichkeits-Transferfunktion zwischen rund r' 9t = genereller Integrationsbereich (1) Wenn man die Quelldichteverteilung anhand der Rauhigkeitslänge einschränkt auf die Höhe der Bodenoberfläche (z = zo) und dem Standort des Messmastes feste Koordinaten zuweist, ergibt sich aus (1) eine Beschreibung der sogenannten Quelldichtefunktion: 4 0000 l1(X m ,Ym ,Zm) = f fQTJ (x',y', z' = Zo)' f(x", -x',y", -y',z", - Zo }dx'dy' (2 ) -00-00 mit: xm, Ym, Zm = Koordinaten des Messmastes x', y', z' = Koordinaten des Quelle In Formel (2) werden alle Messungen auf eine effektive Höhe bezogen, relativ zu einer Verschiebungshöhe, die nicht extra mit aufgeführt wird. Da die Verteilung der Quellstärken nicht im Voraus bekannt ist, muss noch eine weitere Vereinfachung eingeführt werden. Es wird angenommen, dass die Oberfläche aus einer Vielzahl von individuellen Punktquellen besteht, die alle eine einheitliche Stärke haben (unit surface point source). QTJ fließt somit als Konstante QTJ,u in die Berechnungen mit ein. Der unterschiedliche Einfluss jedes dieser einzelnen Punkte auf eine Messung wird nur noch über die Verteilungsfunktion in der Auswertung berücksichtigt. Dadurch lässt sich der relative Anteil jeder Einzelfläche am Gesamteffekt mit dem vertikalen Fluss gleichsetzen, der von dieser ausgeht. Die Verteilungsfunktion fist strenggenommen abhängig von den Eigenschaften des Parameters 11. Deren Transportcharakteristika zwischen Quelle und Sensor unterscheiden sich z.B. weitgehend zwischen radiativ transferierten bzw. turbulent diffundierten Parametern. Betrachtet man aber die Form der Quelldichtefunktion über einer unendlichen Fläche, so ähnelt sie sich für diese beiden Gruppen weitgehend, die Unterschiede können also vernachlässigt werden. Formel (2) kann in der angegebenen Form in der Praxis schlecht eingesetzt werden, da über eine unendliche Fläche integriert werden müsste. Ein besser anwendbares Konzept ist es daher, nicht den gesamten Effekt auf einen gegebenen Standort zu berechnen, sondern eine Region abzuschätzen, die am wahrscheinlichsten den Wert von 11 in der Höhe Zm beeinflusst. Schon die Berechnung einer Fläche, von der 99 Prozent des Bodeneinflusses auf das Messergebnis ausgehen, ermöglicht eine diskrete Einschränkung der source-area über eine festgelegte Außengrenze. Verwendet man diese als Integrationslimit, ist das Integral über die Quelldichtefunktion gleich der Wahrscheinlichkeit, mit der die berechnete Fläche einen Einfluss auf die Messung ausübt. Man kann eine solche Abgrenzung für jede beliebige Effekt-Level P am Gesamteffekt bestimmen. SCHMID & OKE (1990) bezeichnen die kleinste Fläche, die für einen bestimmten Wert von P, beispielsweise für P = 0,5, berechnet werden kann, als Quellgebiet Op für das Niveau P. Diese Fläche wird begrenzt durch die sogenannte Isoplethe der Quelldichtefunktion für P, fp . Der Anteil der gesamten Quelldichtefunktion <Ptob welcher durch Integration innerhalb der Grenzen von fp berechnet werden kann, wird bezeichnet als <pp. P kann demnach über folgende Formel nach SCHMID (1997) berechnet werden: SI f(~,~,/).z)d~d/)'y p = <P p = _-'-p _______ _ <Ptot '" '" f ff(~,/).y,/).z)d~d/).y ( 3 ) mit: /).x = Xm -X = Abstand in Hauptachsenrichtung zwischen Sensor und Quelle /).y = Ym - Y = Abstand quer zur Hauptachsenrichtung zwischen Sensor und Quelle /).z = Zm - Z = vertikaler Abstand zwischen Sensor und Quelle P = Anteil der Quelldichtefunktion innerhalb der Fläche Op am Gesamteffekt <pp = Integral der Quelldichtefunktion innerhalb der Fläche Op <Ptot = gesamte integrierte Quelldichtefunktion Eine räumliche Darstellung eines berechneten Quellgebiets der Bodeneinflüsse auf einen Sensor bietet Abbildung 2. Hier wird auch das Prinzip der Effekt-Level P mit den dieser zugeordneten Parametern Isoplethe <pp und Quellgebiet Op verdeutlicht. Die generelle Form des Quellgebiets bildet durch die im Abbildung 2: Die Quelldichtefunktion und ihre Auswirkungen auf das Quellgebiet der Bodeneinflüsse auf einen Sensor. Nach SCHMID (1994). 5 Vergleich zur lateralen Komponente vorherrschende Verteilung des Parameters in der Windrichtung x eine eiförmige Ellipse. Die Isoplethen <pp für unterschiedliche Effekt-Levels P umschließen mit einem wachsenden Wert von P eine immer größere Quellgebietsfläche üp. Für alle einzelnen source-areas der unterschiedlichen P-Niveaus liegt der Punkt der maximalen Quellstärke fmax im Zentrum. Jedes dieser Quellgebiete kann, da sich ihre generelle Form nicht voneinander unterscheidet, anhand von wenigen charakteristischen Parametern eindeutig beschrieben werden. Diese Kenngrößen werden in der folgenden Abbildung 3 dargestellt. Mit: Xm = Abstand zwischen dem Punkt des größten Einflusses und Sensor gegen die Hauptwindrichtung [km] a = minimaler Abstand zwischen Sensor und Quellgebiet [km] e = maximaler Abstand zwischen Sensor und Quellgebiet [km] d = maximale Breite einer Hälfte des Quellgebiets [km] Abbildung 3: Charakteristische Dimensionen des Quellgebiets der Bodeneinflüsse. Nach SCHMID (1994). Insgesamt kann, wie in Abbildung 3 verdeutlicht wird, das Quellgebiet anhand von 5 charakteristischen Parametern beschrieben werden. Der Wert Xm markiert den Abstand zwischen Sensor und dem Punkt des maximalen Bodeneinflusses auf die Messung. Diese Distanz ist immer entgegen der Windrichtung ausgerichtet. Des weiteren werden der minimale und der maximale Abstand a bzw. e zwischen Sensor und dem äußeren Rand des Quellgebiets zur genauen Definition der Fläche benötigt. Nicht eingezeichnet ist die Entfernung }(,J, welche den Abstand zwischen Sensor und dem Punkt entgegen der Windrichtung markiert, an dem jede Hälfte des Gebiets die maximale Breite d aufweist. Eine letzte wichtige Kenngröße, mit der unterschiedliche Quellgebiete miteinander verglichen werden können, ist deren Fläche, von SCHMID (1994) bezeichnet als A. Dieser Parameter ist vor allem aussagekräftig bei einer Gegenüberstellung der Ergebnisse verschiedener Effekt-Level P. Die für die Berechnung des Quellgebiets notwendigen Eingabeparameter sind die Höhe des Sensors über Grund Zm, die Obuchov-Länge L, die Rauhigkeitslänge 20, die Schubspannungsgeschwindigkeit u. sowie die Standardabweichung der lateralen Windkomponente cr y• Daraus lassen sich mit Zm/zQ , ZmIL und crv/u. drei dimensions lose Skalierungsfaktoren bilden. Obwohl die Rauhigkeitslänge selbst 6 als Eingabeparameter verwendet wird, lässt sich das Verfahren nach SCHMID (1997) auch zu einer footprint-abhängigen zo-Bestimmung verwenden. Dazu ist die Anwendung eines iterativen Verfahrens notwendig, das in den nächsten Abschnitten näher erläutert werden soll. Alle übrigen Eingabegrößen lassen sich mit Hilfe von Parametrisierungsansätzen auf Grundlage einfacher meteorologischer Datensätze abschätzen. 3 Die Bestimmung einer footprint-abhängigen Rauhigkeitslänge Wie im vorhergehenden Abschnitt schon erwähnt steht mit den footprint-Analysen ein geeignetes Instrument zur Verfügung, um den Einfluss einer räumlich differenzierten Geländeoberfläche auf Profilmessungen abschätzen zu können. Von besonderem Interesse bei einer derartigen Untersuchung ist die Rauhigkeitslänge zoo Dieser Parameter ist repräsentativ für die Auswirkungen, die vom Untergrund auf das anströmende Windfeld ausgeübt werden. Bei einem hohen Wert von Zo kommt es aufgrund der Reibung zu verstärkter Turbulenzbildung in den bodennahen Luftschichten. Ideale Standorte für meteorologische Messungen weisen, zumindest in der Hauptwindrichtung, ein ebenes und störungsfreies Gelände auf. Die Bedingungen werden um so optimaler, je geringer die Rauhigkeitslänge auseillt. Gerade im mitteleuropäischen Binnenland, das typischerweise dicht besiedelt ist und häufig die Landnutzungsform ändert zwischen Wäldern und Ackerbzw. Weideflächen, ist es jedoch praktisch unmöglich, ein störungsfreies Anströmungsgebiet zu für hohe Masten zu finden. Weicht man auf weniger optimales Gelände aus, ist es notwendig, den genauen Einfluss der Störfaktoren möglichst schon im Voraus abschätzen zu können. 3.1 Generelles Ablaufschema der Methode In der vorliegenden Untersuchung wird eine Methode vorgestellt, mit der sich die veränderliche Rauhigkeitslänge berechnen lässt, die sich aus dem in Größe und Lage variierenden Quellgebiet der Bodeneinflüsse ergibt. Um die Theorie der footprint-Analysen für einen solchen Anwendungszweck einsetzen zu können, wurde folgende generelle Vorgehensweise entworfen: 1. Da wie im vorhergehenden Abschnitt schon erwähnt die Rauhigkeitslänge einer der Eingabeparameter ist, die für das hier eingesetzte footprint-Modell nach SCHMID (1997) benötigt werden, muss mit einem iterativen Ansatz gearbeitet werden. Die Berechnung startet mit einem Schätzwert für zo, welcher annäherungsweise einen Mittelwert darstellt aus dem Spektrum der zu erwartenden Ergebnisse. In der vorliegenden Untersuchung wurde ein Anfangswert von Zo = 0,1 m gewählt, was nach TROEN & PETERSEN (1989) der allgemeinen Rauhigkeitsklasse 2 für landwirtschaftliche Nutzflächen mit windbeeinflussenden Hindernissen entspricht (s.u.). 2. Um den Einfluss der variierenden Quellgebiete auf eine Windgeschwindigkeitsmessung bestimmen zu können, muss für die Umgebung des Standortes als Berechnungsgrundlage ein Datensatz der Rauhigkeitselemente erstellt werden. Dazu werden die topographischen Karten' des umliegenden Geländes für den Mast gerastert, d.h. ausgehend von einer vom Anwender festzulegenden Kantenlänge LR wird das Gebiet in Einzelflächen LR2 unterteilt. Diesen Quadraten wird jeweils eine feste Rauhigkeitslänge zugeordnet, ihre Koordinaten werden relativ zum Messpunkt definiert, welcher sich im Ursprung des entstandenen Systems befindet. Ein auf diese Weise erstellter Da- * Zur Validierung verwendetes Kartenmaterial: Standort Lindenberg: Topographische Karten 1 : 25000 (A V), Stand 1986, VEB Kombinat Geodäsie und Kartographie, herausgegeben vom Ministerium des Innem, Verwaltung Vermessungsund Kartenwesen 7 tensatz vereinfacht die reell vorliegenden Bedingungen, indem die eigentlich kontinuierlich wechselnden Verhältnisse in eine diskrete Matrix übertragen werden. Zur Bestimmung der flächen gemittelten Rauhigkeitslänge wird ein weiter unten vorgestelltes Mittelungsschema eingesetzt, welches die Abschätzung eines repräsentativen Zo-Werts aus den Karteninformationen ermöglicht. 3. Bevor die footprint-Analyse selbst durchgeführt werden kann, müssen die benötigten Eingabeparameter in Zwischenschritten berechnet werden. Für die Bestimmung der Obuchov-Länge L und der Schubspannungsgeschwindigkeit u. sowie der Standardabweichung der lateralen Windkomponente O"v werden Parametrisierungsverfahren eingesetzt, welche mit Hilfe von einfachen meteorologischen Messungen arbeiten (s.u.). Die Messhöhe Zm, für welche die footprint-Berechnung durchgefuhrt werden soll, kann jeweils frei gewählt werden. Sollen individuelle Rauhigkeitslängen für unterschiedliche Höhenniveaus eines Mastes berechnet werden, erfolgt die Berechnung der Stabilität der atmosphärischen Schichtung lediglich fur die unterste dieser Höhen, da die benötigten Eingabeparameter boden nah gemessen werden. Die dabei ermittelten Ergebnisse werden für die oberhalb liegenden Sensoren übernommen. Die Rauhigkeitslänge wird im ersten Durchgang der Iteration wie oben erwähnt abgeschätzt, bei den weiteren Schritten übernimmt man das Ergebnis des vorhergehenden Durchlaufs. 4. Mit der footprint-Analyse nach SCHMID (1997) kann nun die Ausdehnung des Quellgebiets der Bodenrauhigkeit berechnet werden. Die im weiteren Verlauf dieses Iterationsverfahrens benötigten charakteristischen Lagemaße sind minimaler und maximaler Abstand zwischen Sensor und source-area, a bzw. e, sowie die maximale Breite der Fläche d und deren Abstand zum Sensor entgegen der Windrichtung XJ. Siehe dazu auch Abbildung 3. 5. Das berechnete Quellgebiet wird mit Hilfe der Windrichtungsmessung auf die korrekte Position innerhalb der Rauhigkeitslängen-Matrix projiziert. Anschließend werden die innerhalb dieser Fläche liegenden Rasterpunkte anhand geeigneter Berechnungsalgorithmen identifiziert. Jedem dieser Werte kann mit Hilfe der Quelldichtefunktion eine Gewichtung zugewiesen werden. Die auf diese Weise erhaltenen modifizierten zo-Daten bilden im Mittelwert die effektive Rauhigkeitslänge des vorliegenden Areals. 6. Das in 5. erhaltene, verbesserte Ergebnis von Zo wird in 2. verwendet, um einen neuen Durchlauf der footprint-Analyse zu starten. Diese Iterationsschleife wird so lange wiederholt, bis die absolute Differenz zwischen der in 2. eingesetzten und der unter 5. berechneten Rauhigkeitslänge unter eine bestimmte, vom Anwender festgelegt Schranke fällt. Typischerweise ist die Iteration nach drei Schritten beendet, nur in Ausnahmefalien können auch fünf bis sechs Durchgänge notwendig sein. Die oben beschriebene Vorgehensweise zur Berechnung einer footprint-abhängigen Rauhigkeitslänge wurde in ein FORTRAN-Programm mit dem Titel FOOTSTAB umgesetzt. In den folgenden Abschnitten werden die Methoden und Algorithmen vorgestellt, die zur Erstellung und Anwendung dieses Programms notwendig sind. 3.2 Die Erstellung einer Rauhigkeitslängen-Matrix Informationen über die Geländebedingungen in der Umgebung des Mastes werden wie oben schon erwähnt in eine Matrix der Rauhigkeitslängen umgesetzt. Als Grundlage dazu dienen allein die topographischen Karten der Region, aus welchen nach einem im weiteren Verlauf näher beschriebenen Schema die zo-Daten entnommen werden. In einem ersten Schritt wird der Auflösungsgrad der Rasterung über den Abstand der Punkte zueinander, LR, festgelegt. Jede der abgegrenzten quadratischen 14 q sal = 0,622. e sal (16 ) p e =: 6112. ex [17,67, (T -273,16)] sal' P T - 29 66 , ( 17) mit esat = Sättigungsdampfdruck [hPa] p Luftdruck [hPa] Die Formeln (16) und (17) nach STULL (1988) können auch verwendet werden, um die beiden Feuchte-Parameter des Faktors Fw aus Gleichung (12) zu bestimmen. Für die spezifische Feuchte bei Sättigung qsat kann genauso verfahren werden wie oben angegeben, zur Berechnung der spezifischen Feuchte der Luft qair verwendet man in (16) den Messwert des aktuellen Dampfdrucks anstelle des Sättigungsdampfdrucks. Der turbulente, atmosphärische Widerstand ra im Nenner von (12) ist eine Funktion der Stabilität der Schichtung. Er lässt sich berechnen durch eine Umstellung der Formel (7) unter Einbeziehung der Verschiebungshöhe d: ( Z-dJ In -- z~ - -\jf M ra = (18 ) K ·u. mit: d = Verschiebungshöhe [m] Der Verdunstungswiderstand der Pflanze rp ist ein Parameter, dessen exakte Bestimmung die Kenntnis zusätzlicher pflanzenphysiologischer Parameter erfordert. Eine mögliche Berechnungsformel nach GARA TT (1992) schlägt folgende Abschätzung vor: r =: P LAI mit: rsi = Stomata-Widerstand der Einzelblätter [s/m] LAI:::: Leaf-Area-Index [-] ( 19 ) Beide Parameter können bei genauer Kenntnis von Art, Zusammensetzung und Alter der Vegetation aus Tabellen entnommen werden. Für den Blattflächenindex besteht zudem die Möglichkeit einer direkten Messung bzw. der Abschätzung aus Beobachtungen. Da derartige Angaben im Normalfall nicht zur Verftigung stehen, kann der Verdunstungswiderstand auch als Konstante mit in die Berechnung eingehen. Ein geeigneter Mittelwert ist nach JACOBSON (1999) 60 s/m. 3.3.3 Mathematische Umsetzung Die Anwendung dieser Formelsätze zur Stabilitätsbestimmung, also die Berechnung turbulenter Flüsse aus Profilmessungen mit Hilfe der universellen Funktionen, ist mathematisch anspruchsvoll. Will man beispielsweise die Schubspannungsgeschwindigkeit u. mit Hilfe von Formel (7) ermitteln, benötigt man auf der rechten Seite der Gleichung fur die Bestimmung der integrierten universellen Funktion die Obuchov-Länge L. Diese wiederum hängt von u. ab, so dass ein Zirkelschluss entsteht, der nur iterativ zu lösen ist. In diesen Iterationsprozess mit eingezogen werden muss zusätzlich noch die Berechnung des ftihlbaren Wärmestroms (s.o.). Eine Umsetzung dieser Problematik in ein einfaches, iteratives Berechnungsverfahren entwickelte HOLTSLAG (1984). Der Autor beginnt seine Berechnungsschleife mit einer abgeschätzten Obuchov- 15 Länge. Sein Vorschlag ist ein Startwert L = 00. Eigene Erfahrungen mit dem Verfahren zeigen aber, dass die Algorithmen stabiler arbeiten, wenn man im labilen Fall mit L = -1000 beginnt, im stabilen mit L = 1000. Anhand dieses Schätzwertes kann durch Umstellung der F1uss-Profil-Beziehung (7) eine erste Schubspannungsgeschwindigkeit berechnet werden: (20 ) Unter Verwendung Schubspannungsgeschwindigkeit aus (20) und des turbulenten fühlbaren Wärmestroms aus (11) ergibt sich aus Formel (5) ein neuer Wert für die Obuchov-Länge L. Dieser wird wiederum eingesetzt in die Formeln (11) und (20), um verbesserte Werte der Skalierungsfaktoren zu erhalten. Der oben beschriebene Iterationsprozess wird so lange fortgesetzt, bis sich der eingesetzte Wert der Obuchov-Länge dem Berechnungswert auf die geforderte Genauigkeit angenähert hat. Laut HOLTSLAG (1984) sind im allgemeinen nicht mehr als 3 Schritte notwendig, um die Fehlergrenze von 5 Prozent zu unterschreiten. 3.3.4 Benötigter Datensatz Unter Einbeziehung der Parametrisierung der Strahlungsbilanz aus Messungen der Globalstrahlung benötigt das beschriebene Verfahren somit unter Einbeziehung des Penman-Monteith-Modells folgenden meteorologischen Größen zur Berechnung des fühlbaren Wärmestroms: Globalstrahlung [Km/s] Lufttemperatur [K] Dampfdruck [hPa] Luftdruck [hPa] Windgeschwindigkeit [m/s] Zusätzlich sind noch einige Parameter zur Standortbeschreibung notwendig: Geographische Position des Standorts Rauhigkeitslänge und Verschiebungshöhe [m] Angabe von Uhrzeit und Datum 3.4 Die Standardabweichung der lateralen Windkomponente Die Standardabweichung der lateralen Windkomponente O"v kann mit Hilfe des BodenschichtSkalierungsparameters zlL abgeschätzt werden über folgende Formel nach FOKEN et al. (1997): 0"1'=2,7·u. ,fürl~I<0,03l9 (21) ", ~ 4,15 ·u •. (~t ,flir 0,0319 < I~I< 1 (22) mit: O"y = Standardabweichung der lateralen Windkomponente [m/s] Die Beziehungen (21) und (22) wurden ursprünglich entwickelt, um die Standardabweichung der Strömung O"u in Windrichtung abschätzen zu können. Nach WYNGAARD & CLIFFORD (1978) ergeben sich aber die gleichen Zusammenhänge auch für die Streuung der lateralen Windkomponente, 16 so dass die Formeln in der oben aufgefUhrten Darstellung für den angestrebten Verwendungszweck einsetzbar sind. 3.5 Die Kombination von source-area und Rauhigkeitslängen-Matrix Mit den im vorhergehenden Abschnitt aufgefUhrten Formelsätzen können alle Eingabeparameter fUr das footprint-Modell nach SCHMID (1997) berechnet werden. Mit Hilfe dieser Ergebnisse kann fUr eine festgelegte Sensorhöhe Zm das Quellgebiet der Bodenrauhigkeit mit dem Verfahren nach SCHMID (1997) bestimmt werden. Als Ausgabeparameter dieses Moduls des Programms FOOTSTAB werden die charakteristischen Lagemaße der source-area a, e, d und Xd sowie die Lage des Punktes des maximalen Einflusses, Xmax, ausgegeben. Im diesem Abschnitt soll nun ein Algorithmus vorgestellt werden, mit dem es möglich ist, das anhand dieser Lagemaße bestimmte Quellgebiet mit der zo-Matrix zu verschneiden. Ziel des Verfahrens ist es, fUr eine festgelegte Fläche die innen liegenden Punkte des vorliegenden Rasters zu identifizieren. Aus diesen mit Hilfe der Quelldichtefunktion ein gewichteter Mittelwert der effektiven Rauhigkeitslänge bestimmt. Um die Matrixpunkte innerhalb des Quellgebiets bestimmen zu können, muss zuerst eine geeignete Formel gefunden werden, die dessen äußeren Rand nachvollziehen kann. Als weitere Annahme werden die parallel zur Ordinate übereinander liegenden Punkte des zo-Rasters für die Berechnung als durchgehende Gerade behandelt. Das Grundprinzip der folgenden Algorithmen ist es nun, Schnittpunkte zu identifizieren zwischen der Grenze des Quellgebiets und diesen einzelnen Geraden der Rauhigkeitsmatrix. Dabei können 2 Fälle auftreten: Wird keine Überschneidung von Gerade und der Formel fUr den Rand der source-area gefunden, so spielt die betrachtete senkrechte Reihe von Matrixpunkten keine Rolle für die weiteren Berechnungen. Ergeben sich jedoch 2 Schnittpunkte (ein einzelner Schnittpunkt ist nicht möglich), so bildet die Gerade eine Sekante des Quellgebiets der Bodeneinflüsse. In diesem Fall muss weiter überprüft werden, ob zwischen den beiden gefundenen Koordinatenpunkten Matrixpunkte auf der betrachteten Linie und damit auch innerhalb der source-area liegen. 3.5.1 Herleitung der Formel für die Quellgebiets-Grenzen Wie weiter oben schon erwähnt ergibt sich aus der footprint-Analyse nach Schmid ein Quellgebiet, dessen Form generell als eiförmig bezeichnet werden kann. Die Formel fUr deren äußeren Rand orientiert sich daher grundsätzlich am Einheitskreis. Dieser berechnet sich nach folgender Grundgleichung: 2 2 1 Xo + Yo = (23 ) mit: XQ, Yo = Koordinaten des Einheitskreises Da das Quellgebiet vor bzw. jenseits des Punktes der maximalen Breite, Xd, vom Maststandort aus betrachtet eine unterschiedliche Form ausweist, muss die Fläche fUr die Berechnung in 2 Hälften unterteilt werden. Im weiteren Verlauf der Herleitung einer Formel fUr den äußeren Rand der source-area wird daher differenziert zwischen einer sensomahen und einer sensorfemen Ellipsenhälfte. Die Windrichtung wird vorerst nicht berücksichtigt, d.h. die Ausrichtung der Fläche erfolgt entlang der positiven Ordinate, entsprechend einem Wind direkt aus Norden. Für die sensomahe Ellipsenhälfte ergibt sich, aufgelöst nach y, folgende Grundformei: Yo = -Jlx~ (24 ) Die sensorfeme Ellipsenhälfte berechnet sich ausgehend von folgender Grundgleichung: 17 Yo := Jl-xg (25 ) Im nächsten Schritt müssen die vom footprint-Modell ausgegebenen charakteristischen Lagemaße a, e, d und X<J dazu verwendet werden, aus den beiden bislang noch runden Kreishälften jene Halb-Ellipsen zu formen, die dem mit dem footprint-Modell nach SCHMID (1997) berechneten Quellgebiet entsprechen. Dazu werden rur die sensorferne Ellipsenhälfte folgende Beziehungen verwendet: x xo·d (26) Y=Yo ,{x d -a)+x d (27) mit: d maximale Breite des Quellgebiets [km] X<i Entfernung des Punktes der maximalen Breite vom Sensor in Windrichtung [km] a minimale Entfernung des Quellgebiets vom Sensor [km] Nach Einbeziehung dieser Parameter in Formel (24) ergibt sich rur die sensornahe Ellipsenhälfte folgende Beziehung: r)_(~)2 .(x, -al +x, (28 ) Die Berechnungen ftir die sensorferne Ellipsenhälfte haben eine leicht abweichende Form: x xo·d (29) Y=Yo.(e-xd)+x d (30) mit: e == maximale Entfernung des Quellgebiets vom Sensor [km] Die Einbeziehung von (29) und (30) in (25) ergibt folgende Formel ftir die sensorferne Ellipsenhälfte: (31 ) Eine schematische Darstellung der oben beschriebenen Vorgehensweise bietet Abbildung 6. 3.5.2 Rotation des berechneten Quellgebiets in die Windrichtung Durch eine Vektorisierung der Umfangskoordinaten dieses Quellgebiets ist es möglich, die umrandete Fläche entsprechend der Windrichtung um den Ursprung rotieren zu lassen. Auf diese Weise kann die Lage des eigentlichen Quellgebiets berechnet werden, denn bei den bisherigen Berechnungen war die Windrichtung auf Nord festgelegt. Der erste Schritt dazu ist die Berechnung des Abstandes eines Punktes vom Ursprung, also der Länge des Ortsvektors. VL ~X2+ y2 (32 ) mit: VL = Vektorlänge Die Berechnung des Winkels dieses Vektors gegen Nord erfolgt über: a ~arcWn(; J (33 ) mit: a = Abweichung des Längenvektors eines Punktes gegen Nord [rad] 18 Obere Hälfte des Einheitskreises untere Hälfte des Einheitskreises ........ ............. ........................... 3 ., ..................... , L ........................... ... L.............. ..., ................ .:3 . .. L ............................ L .................... . ... . .. L .. . ...... . ..... • ....... ............................ ..... ....................... . .:3 .J .......... .. .. ... L ...... ... ........... ... ... . L .... ..................... .1 Einbeziehung der charakteristischen Dimensionen des Quellgebiets: a, Xd, e, d (hier: a = 1; Xd = 4; e = 9; d = 2) ............ ....................................................................... 10 .,. ······· r·············.,. ............ , 10 . R / \ 7 I \ . \. , '- -5 .. .:3 ·2 ·1 o 1 2 5 -5 ~ -3 ~ zusammengesetzt .... ... '0 . , ...................... . ~ --- / / n -, _R _R ................................................... c............ .. ............ . .. -10 ... J ./ Abbildung 6: Schematische Darstellung zur Herleitung der Funktion zur Beschreibung des Quellgebiet-Rands. Die Bestimmung der neuen, im die Windrichtung rotierten Koordinaten erfolgt unter Einbeziehung von (32) und (!Textmarke nicht definiert, NEU7) über folgende Beziehungen: X rot = VL . sin(a +~) Yrot = VL . cos(a +~) mit: Xrot = x-Koordinate des um die Windrichtung ~ rotierten Punktes Yrot = y-Koordinate des um die Windrichtung ~ rotierten Punktes ~ = Windrichtung [rad] (34 ) (35 ) 19 Eine mit Hilfe der Formeln (32) bis (35) durchgeführte Rotation des berechneten Quellgebiets (Abbildung 6) ist in Abbildung 7 dargestellt. Die Windrichtung beträgt hier 1680 • . . ..................... ................... _ .. . 1... ....................... .1.. .......................... " . .• ...................... 1.. .... .... ................... 1 Abbildung 7: Quellgebiet aus Abbildung !Textmarke nicht definiert, FOOT6, rotiert um 1680 Eine Rotation des berechneten Quellgebiets auf das tatsächlich interessierende Areal in Windrichtung ist demnach mit einfachen Formelsätzen durchführbar. Die mögliche Verschneidung mit der erstellten Rauhigkeitslängen-Matrix ist in der folgenden Abbildung 8 dargestellt. Auch hier ist das Quellgebiet, nach der Transformation des Einheitskreises in die Ellipsenform, in eine angenommene Windrichtung von 168 0 gedreht worden. Wie weiter oben schon erwähnt werden jeweils parallel zur Ordinate übereinander liegende Rasterpunkte, hier als Kreuze (x) dargestellt, als eine durchgehende Gerade aufgefasst. Im Bereich des Quellgebiets bildet diese Gerade eine Sekante. Sind die beiden Schnittpunkte bekannt, kann anhand eines entworfenen Berechnungsschemas festgestellt werden, welche der Rasterpunkte auf dieser Geraden innerhalb des Quellgebiets liegen. Hervorgehoben ist in !Textmarke nicht definiert, FOOTIA6 die Linie durch einen Abszissenwert von x= 1. x x x x x .- 1 :ek an le ~ urch x = 1 I : / I ~ Ch"' ft Pun kt 1 I: Abbildung 8: Darstellung der Schnittpunktberechnung zwischen rotiertem Quellgebiet der Bodenrauhigkeit und einer ordinaten parallelen Sekante durch die Abszisse bei x = 1. Die in diesem Abschnitt vorgestellte Berechnungsmethode ist allerdings weniger geeignet fur den eigentlichen Zweck der Untersuchung, die Identifizierung der innerhalb des Quellgebiets liegenden Rasterpunkte der Rauhigkeitslängen-Matrix. Das grundlegende Konzept für einen geeigneteren Algorithmus wird im folgenden Abschnitt erläutert. 3.5.3 Ein Algorithmus zur Verschneidung von Quellgebiet und zo-Matrix Die in Abbildungen 7 und 8 dargestellte Rotation des Quellgebiets selbst ist, wie oben schon erwähnt, nicht die optimale Variante zur Schnittpunktbestimmung zwischen deren äußerem Rand und den ordinatenparallelen Geraden der zo-Matrix. Die mathematische Umsetzung der Identifikation von inner- 20 halb des Quellgebiets liegenden Rasterpunkten kann leichter durchgeführt werden, wenn man die nach Norden orientierte ellipsenförmige Fläche aus Abbildung 6 verwendet, und in diese die um den Ursprung rotierten Sekanten legt. Die Steigung der rotierten Sekanten ist somit abhängig von der Windrichtung. Zusätzlich muss noch eine Verschiebung entlang der Ordinate vorgenommen werden, deren Wert von der relativen Position der Sekante im Koordinatensystem (s.u.) berechnet wird. Diese Anpassung kann über die folgende Geradengleichung erreicht werden: x S·L Y= + __ R tan(7t -~) sin~ mit: S = laufende Nummer der Sekanten (Ursprungsgerade = 0) LR = Abstand der Rasterpunkte voneinander [km] (36 ) Der Verlauf zweier rotierter Sekanten durch das oben schon vorgestellte Quellgebiets-Beispiel (Windrichtung = 168°) ist in der folgenden Abbildung 10 dargestellt. V . IO I ~ö.a;;e~;y .... .. ------ / / /y\: / ekan e S = 01 . 11 ISteigung = 1Itan(l"C ) 1 17 II AbS!~=us,"~ \ I~ 1/ V / / '----- r--. I / --- V /n / -3 -2 -1 Abbildung 9: Schematische Darstellung zweier in das Quellgebiet hineinrotierter Sekanten. Zur Berechnung der Schnittpunkte zwischen den einzelnen Sekanten und dem Rand des Quellgebiets müssen die Formeln (28) und (36) für die sensornahe Ellipsenhälfte bzw. (31) und (36) für die sensorferne Ellipsenhälfte gleichgesetzt werden. Für die Abszissenpunkte des sensornahen Anteils der Ellipse ergibt sich daraus durch die Auflösung nach x: K~ 1 2 X(U)12 =± ---+E u -Eu , Fu Fu (37 ) Die Komponenten Ku, Fu und Eu ergeben sich dabei als Zwischenergebnisse aus den charakteristischen Lagemaßen des Quellgebiets sowie den oben schon eingeführten Parametern Sund LR, welche die Eigenschaften der Rauhigkeitslängen-Matrix näher definieren. (38 ) (39 ) ( 40) 21 Für die sensorfeme Ellipsenhälfte gelten folgende Beziehungen: ( 41 ) Die Eingabeparameter 1<0, Fo und Eo berechnen sich nach ähnlichen Formeln wie die oben in (38) bis (40) schon aufgestellten Beziehungen für den sensomahe Anteil des Quellgebiets: K _ S· LR xd 0sin ~ . (e -x d ) (e -x d ) ( 42) Fo =-(~)' {an(" -~~.(e-XJJ ( 43 ) (44 ) Die genaue Herleitung der Formeln (37) bis (44) wird im Anhang beschrieben. Die y-Koordinaten der Schnittpunkte lassen sich über Gleichung (36) bestimmen. Falls der Wurzelausdruck in (37) bzw. (41) negative Werte annimmt, existieren keine Schnittpunkte zwischen der betreffenden ordinatenparallelen Geraden der zo-Matrix und dem äußeren Rand des Quellgebiets. Im FOOTSTAB-Programm werden in diesem Fall die Ergebnisse gleich Null gesetzt, und die Berechnung wird abgebrochen. Ansonsten ergeben sich aus diesen beiden Formeln insgesamt 4 Schnittpunkte für die gesamte Ellipse: x(u)], X(U)2, x(o)] und X(O)2' Daher müssen die auf diese Weise ermittelten Ergebnisse noch auf ihre Plausibilität untersucht werden, da reell nur 2 Schnittpunkte vorhanden sind. Grundsätzlich müssen dabei die Ordinatenwerte der berechneten Schnittpunkte innerhalb der durch die charakteristischen Lagemaße vorgegebenen Intervalle des jeweiligen Abschnitts des Quellgebiets liegen. Dabei gilt für die sensomahe Ellipsenhälfte a<y(u)],2<Xd, für die sensorfeme Xd<y(o)],2<e. Auf diese Weise lassen sich mit den vom footprint-Algorithmus ausgegebenen Längen a, ~ und e die gesuchten 2 Schnittpunkte eindeutig identifizieren. Die Berechnung der Koordinaten der rotierten Schnittpunkte kann wiederum mit Hilfe einer Vektorisierung vorgenommen werden. Die Vektorlänge wird dabei über Formel (32) ermittelt, der Winkel ergibt sich aus: a ro' = arctan(; ) +~ ( 45 ) mit: arat = Winkel des Schnittpunkt-Vektors gegen Nord im rotierten Koordinatensystem [rad] Die eigentlichen Schnittpunkte des Randes des Quellgebiets mit der jeweiligen Sekante berechnen sich aus diesen Angaben über folgende Formel: SP = VL • cos(a rOI) ( 46) mit: SP = Ordinatenwert des Schnittpunkts im rotierten System [km] Auf diese Weise ergeben sich für jede Sekante 2 Schnittpunkte, die im folgenden Abschnitt dazu verwendet werden, die innerhalb des Quellgebiets liegenden Rasterpunkte zu identifizieren und somit eine effektive Rauhigkeitslänge zu ermitteln. 22 3.5.4 Berechnung einer effektiven Rauhigkeitslänge 3.5.4.1 Identifizierung der innerhalb des Quellgebiets liegenden Rasterpunkte Das Ergebnis der oben beschriebenen Schnittpunkt-Berechnung ist im FOOTSTAB-Programm eine Matrix, die für jede Spalte der Rauhigkeitslängen-Matrix einen oberen und einen unteren Schnittpunkt enthält. Anband dieser Schnittpunkte wird im folgenden überprüft, welche der Matrixpunkte innerhalb des Quellgebiets liegen, indem folgende Formel durchlaufen wird: SP I 5, M(zo) ' . 5, SP 2 I , I ,} I , mit: i = Laufindex der Spalten der Matrix j = Laufindex der Zeilen der Matrix SPi ,1 = oberer Schnittpunkt der i-ten Sekante SPi,2 = unterer Schnittpunkt der i-ten Sekante M(zo)ij = Matrixpunkt an der Position i, j (47 ) In Abbildung 12 wird das zurück in seine Ausgangslage rotierte System von Quellgebiet und Sekanten sowie die vorgenommene Differenzierung der Rasterpunkte dargestellt. ........... >< .. -- ....... ,)( ... · ··· ····· x ........... .. )( ............ ,)( .. . ........ .. Jo:< X )( x)( x x x * ElElElElElEJEJ~ i x * X: :i i i: i ~ x x x * x x )( )( x X .. - ...... ·x· .... · ~tO .. .. .... )C . . ...... )( ............ ~ .......... )( ............ )( Abbildung 10: Schematische Darstellung des in die Windrichtung zurückrotierten Systems von Quellgebiet und Sekanten. Verdeutlichung der Differenzierung zwischen Rasterpunkten innerhalb (0) und außerhalb (x) der source-area. Sind alle Rasterpunkte, die innerhalb der betrachteten Fläche liegen, identifiziert worden, werden für insgesamt 10 verschiedene Effekt-Level des Quellgebiets alle zo-Werte aufsummiert und die Gesamtanzahl der Punkte gespeichert. Der ungewichtete Mittelwert fur die gesamte source-area kann berechnet werden, indem man von der 100%-Effekt-Level diese Gesamtsumme von Zo dividiert durch die Summe der Rasterpunkte. 3.5.4.2 Berechnung einer effektiven Rauhigkeitslänge unter Einbeziehung der Quelldichtefunktion Um auch die Quelldichtefunktion mit in die Berechnung der effektiven Rauhigkeitslänge einfließen lassen zu können, werden, wie oben schon erläutert, die Werte der unterschiedlichen Effekt-Level P verwendet. Siehe dazu auch Abbildung 2. Das footprint-Modell nach SCHMID (1997) kann grundsätzlich für jedes P-Niveau eine individuelle Fläche bestimmen, die mit der gewählten Wahrscheinlichkeit einen Einfluss auf die Messung hat. Standardmäßig gibt dessen FORTRAN-Routine, welche als Modul mit in FOOTSTAB integriert wurde, die Effekt-Level von 10 bis 100 Prozent des Gesamteffekts in Intervallabständen von 10 Prozent aus. Insgesamt werden also fur 10 verschiedene Niveaus 23 unterschiedliche Quellgebiete der Bodeneinflüsse berechnet. Bezieht man die Einzelflächen jetzt nicht auf die Wahrscheinlichkeit, mit der sie den Sensor beeinflussen, sondern auf ihren Anteil an der Auswirkung des gesamten Quellgebiets auf die Messung, ergibt sich folgendes Bild: Die Rasterpunkte, die sich beispielsweise innerhalb der 10-Prozent-source-area befinden, tragen insgesamt ein Zehntel zum Gesamteffekt mit bei. Diese Fläche wird vergleichsweise klein sein, und somit ist auch die Anzahl der innen liegenden Punkte normalerweise stark eingeschränkt. Die Ellipse zieht sich eng um den Punkt des maximalen Einflusses auf den Sensor. Das jeweils nächstgrößere Quellgebiet, von welchem hier 20 Prozent des Gesamteffekts auf den Sensor ausgehen, schließt die gesamten Flächen der kleineren Effekt-Level (in diesem Fall nur die 10-Prozent-source-area) mit ein, erfasst aber normalerweise zusätzlich noch weitere Rasterpunkte. Entfernt man nun aus der Gesamtmenge der fur das 20-Prozent-Niveau gefundenen Punkte jene der 10-Prozent-Fläche, so verbleibt räumlich gesehen ein Ring, dessen Anteil am Gesamteffekt nur noch 10 Prozent ausmacht. Als Berechnungsgrundlage für die gewichtete Mittelung der innerhalb des gesamten Quellgebiets liegenden Rasterpunkte werden die beiden am Ende des vorhergehenden Abschnitts eingeführten Parameter, Gesamtsumme der Rauhigkeitslänge für ein Quellgebiet sowie die Anzahl der innerhalb liegenden Rasterpunkte, für jede der 10 Effekt-Level-Flächen zusammengestellt. Bei stark labiler Schichtungslage kann es vorkommen, dass die Gebiete mit niedrigem P-Niveau so klein sind, dass sie keinen Matrix-Punkt mit einschließen. Um sicherzustellen, dass keine Einzelfläche die gemittelte Rauhigkeitslänge Null zugewiesen bekommt, wird in diesem Fall ein Korrekturverfahren eingesetzt. Dieses ordnet der Fläche jenes Rasterelement zu, welches den geringsten räumlichen Abstand zur Außengrenze aufweist. Anschließend werden aus jedem einzelnen Quellgebiet die Punkte entfernt, die auch innerhalb der nächstkleineren source-area liegen. Dazu werden lediglich Punktanzahl bzw. Rauhigkeitslängen-Summe voneinander subtrahiert. Nach diesem Schritt verbleiben das ursprüngliche 10Prozent-Quellgebiet sowie 9 konzentrische ellipsenförmige Ringe, die dieses in immer größer werdenden Abständen zur maximum-source-location umschließen. Dabei wird die Anzahl der innerhalb jeder Einzelfläche liegenden Punkte mit steigendem Effekt-Level-Wert zunehmen. Durch Division der über die Subtraktion modifizierten zo-Summe durch die verbleibende Anzahl der innen liegenden Rasterpunkte fur jedes Areal ergeben sich 10 Einzelwerte der gemittelten Rauhigkeitslänge. Da jede der so voneinander abgegrenzten Flächen nun unabhängig von seiner Größe ein Zehntel des Gesamteffekts auf die Messung beinhaltet, kann aus diesen 10 Zwischenergebnissen durch einfache arithmetische Mittelung die gesamte Rauhigkeitslänge unter Einbeziehung der Quelldichtefunktion ermittelt werden. Wie unter Abschnitt 3.1 schon erwähnt, wird nach der Berechnung dieser gewichteten effektiven Rauhigkeitslänge geprüft, inwiefern sich dieses Ergebnis vom zo-Eingabewert fur die Berechnungen unterscheidet. Vom Anwender muss an dieser Stelle ein Schwellenwert festgelegt werden. Unterschreitet die Differenz die angegebene Schranke, so kann das Verfahren abgebrochen werden, der ermittelte Rauhigkeitswert stellt das Endergebnis dar. Ist der Unterschied zwischen Startwert und Berechnungsergebnis zu groß, wird die Iteration fortgesetzt. Schematisch ist der Ablauf der FOOTST ABAlgorithmen in Abbildung 13 aufgeführt. 4 Bewertung des erstellten Berechnungsschemas Mit dem FOOTSTAB-Programm konnte ein zuveflässiges Berechnungsverfahren erstellt werden, das in der Lage ist, auf einfacher Datengrundlage über die Bedingungen der Unterlage eine effektive Rauhigkeitslänge zu bestimmen. Als Eingabegrößen werden neben der zo-Matrix lediglich die oben schon aufgeführten Parameter zur Bestimmung des fuhlbaren Wärmestroms aus der Globalstrahlung sowie