Abstract Since the early sixties continuing into the new millenium Egypt has settled on a strategy of reclaiming dessert land through groundwater abstraction, in order to increase their agricultural production. What effects, this strategy has on the Nubian Aquifer underlying large parts of the Western Dessert, has hitherto been unknown. The Gracity Recovery And Climate Experiment (GRACE) twinsatellite mission might be a means to examine this effect, since the mission is designed to measure changes in groundwater in aquifers on a scale of a few hundred kilometers and up. A total of 102 temporally unevenly spaced GRACE datasets covering the period from 2002 to 2010 were aquired from Center for Space Research, University of Texas consisting of a set of Stokes coefficients to degree and order 100. These were then gaussian filtered using a halfwidth of 0km, and averaged over a region of ∼300×300km. The time-series show a seasonal signal, but no interannual trend. This could indicate, that a drawdown is taking place in the summer as a result of irrigation, whereas groundwater is running into the area in the winther thus recharging the aquifer.
Acknowledgements The author of this report would like to express his gratitude to the following people, whom without their support this project wouldn’t have been what it is today. First of all thanks to my supervisor Niels Schrøder for support and guidance through the writing process and for making the field trip to Egypt in January 2011 possible. Also thanks to senior researcher at DTU Space Ole B. Andersen for his support in the implementation of the data processing. The possibillity to discuss data processing issues with a more experienced researcher has been most valuable. Lastly thanks to Dr. Mohsen A. Gameh Ali from University of Assiut for organizing the field trip to Egypt, which was a great learning experience. Aslo the continous support with respect to geographical information about the New Valley area in the following writing process was most valuable.
i
Contents 1 Introduction 1.0.1 Problem formulation . . . . . . . . . . . . . . . . . . . . . 1.0.2 Target group . . . . . . . . . . . . . . . . . . . . . . . . .
1 3 3
2 Study area 2.1 Geology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Climate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 4 10
3 Method
12
4 Theory 20 4.1 Technical description of the GRACE mission . . . . . . . . . . . 20 4.2 Introduction to potential theory . . . . . . . . . . . . . . . . . . . 22 4.3 The disturbing potential . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 Relationship between the stokes coefficients and the earth surface density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 4.5 Handling of GRACE error . . . . . . . . . . . . . . . . . . . . . . 26 5 Results and conclusion
29
Bibliography
31
A Implementation in MatLab
34
ii
Chapter 1
Introduction ”The only matter that could take Egypt to war again is water.” (Former Egyptian president Anwar Sadat, 1979) Located in the center of the Middle East Egypt plays a key role in the development of the region. Especially since the middle east is known as a politically unstable region, the situation in Egypt is worth studying. Egypts has fought two wars against Israel, and other countries in the region are notoriously known to work for a politically unstable climate. Moreover, water scarcity, lack of development, and too low food production are problems influencing the whole area. Because of this the region has attracted much interest from international development organisations and international agencies during the last 20 years or more. Egypt is located in one of the driest regions of the world (Popp, 2004, p. 17). The Nile flowing from Sudan and Ethiopia is the primary source of water, and by far the largest part of the argicultural production is taking place in the Nile Valley. The Nile Valley is only occupying 4% of the countrys area, the rest being desert (Ibrahim and Ibrahim, 2003, p. 30). The amount of available water for agriculture in Egypt is regulated by the Nile agreement from 1959, assigning 55.5 billion cubic meters of water a year to Egypt (Ibrahim and Ibrahim, 2003, p. 75). Egypt is as well located in a climate with very little rain. The average annual rainfall is 200mm at the Mediterranean coast declining to approximately 0mm south of Cairo (Amer et al., 2005, p. 1). This means, that almost all the agriculture is irrigated. From the point of view of the average egyptian, the rising food prices combined with the widespread poverty is a main problem. The poverty is a result of a lack of industrialization despite many donor programs attempts to counteract the situation (Mitchell, 1995). Also the large corruption plays a role in these problems preventing development and good governance (Ibrahim and Ibrahim, 2003, p. 7). Moreover, the unequal distribution of land means that many egyptian farmers are living in great poverty. They can’t afford to buy food due to the high food prices on the international market, where subsidized grains from EU and USA are being sold, moreover, their land areas aren’t large enough to support the farmers and their families due to the unequal distribution of land. This is also related to the fact, that Egypt during the 90’s pursued an agricultural strategy switching from grain production to animal production due to a 1
growing middle class in the cities thus importing especially wheat. Moreover, it is difficult for the country to access for instance the market in EU due to custum barriers and international trade conventions (Mitchell, 1995). Seen from the governments point of view the agricultural production is too low. This is so, because Egypt has to import massive amounts of primarily wheat every year leading to large state deficits and a large foreign debt (Ibrahim and Ibrahim, 2003, p. 128). As a means to solve this problem, former president Mubarak decided in 1997 to pursue the strategy initiated by president Nasser about reclaiming a part of the dessert to use for agricultural land (Meyer, 2004, p. 195). This is also claimed to be part of the solution to the huge unemployment in the country (Adriansen, 2009). Again this is connected to corruption since the initiative is lower in a corrupt society. Since the rule of president Nasser, Egypt has had a large public sector. Nasser believed in socialism, and a large state should be pulling the development of the country. Unfortunately this has lead to inefficiency in the state apparatus since too many people are employed with little work to do. Due to the lack of social security in Egypt it is not possible to lay off these superfluous employees (Ibrahim and Ibrahim, 2003, p. 93). This is also related to the lack of democracy in the country meaning that laying off public employees could lead to riots. The lack of social security also means that the birth rate is large. This is due to the fact, that the daugthers go to live with their husbands family when they get married. When the parents get old, they have to have some sons to take care of them and since the infant mortality is high, the birth rate also becomes high in order to counteract that some of the children may die (Mitchell, 1995). One way to assess the land reclamation strategy could be to analyze the sustainabillity of the strategy. This becomes especially relevant since the national water resources policy from 2002, with a planning horizon from 1997 to 2017, states that the management of the water is supposed to be ”sustainable” from a socio-economic and environmental point of view (Amer et al., 2005, p. 13). One way to define sustainabillity could be the definition from the Brundtland report: ”Development that meets the needs of present generations without compromising the abillity of future generations to meet their own needs” (World Commission on Environment and Development, 1987). Taking point of departure in the concept of weak sustainabillity, a sustainabillity assessment would consist of an assessment of the natural capital, the man-made capital, and the human capital. A prognosis for the development of these over time would be set up and through an adequate choice of pricing mechanism, it would be possible to calculate whether the development is sustainable. Based on this measure of sustainabillity political choices of meaures to make the development more sustainable could be made. In this way, assessing the sustainabillity of the land reclamation strategy, could serve as a basis for planning and regulation. The natural capital in the Western Dessert is strongly connected to the amount of available water and part of assessing the sustainabillity of the strategy will therefore be assessing the amount of water. The amount of water in the future is dependent on the recharge and discharge characteristics of the studied area. Ahmed (1999) reviews various theories for the recharge of the aquifer under consideration. Some theories argue that the aquifer is completely sealed and that the current groundwater pumping is mining, whereas other theories argue for a certain amount of recharge. It is thus the aim of this project to contribute to a resolution of this discussion. 2
The Western Dessert of Egypt is a large area with very few measurement stations. Information about the water storage is thus in this project inferred from the GRACE twin-satellite. Since little research on water management based on GRACE-data has been performed assessing the applicabillity of using GRACEdata as a tool in water management also becomes the aim of this report. Lastly the aim of the report is to be an introduction to the study of groundwater from space, thus this report will strive to explain most concepts and calculations in detail.
1.0.1
Problem formulation
This has lead to the following problem formulation: What is the magnitude of the temporal development in the water storage of the aquifer underlying Kharga oasis and Dakhla oasis1 in the Western Dessert of Egypt?
1.0.2
Target group
The target group of this report is expected to have a fundamental knowledge of satellite remote sensing. Moreover, readers are also expected to have a minimum knowledge of groundwater theory. Chapter 4 draws on andvanced mathematics, and the details of this chapter will prerequisite knowledge of spherical coordinates, vector calculus, and imaginary numbers. In this respect knowledge of classical mechanics will also be beneficial. The target group of this report will primarily be researchers working in the field of groundwater and/or remote sensing. Since this report collects the theory behind the GRACE mission, it will function as a stepping stone into studying groundwater from space.
1 See
section 2 for a map of the study area.
3
Chapter 2
Study area In this chapter, the current situation in the area with respect to geology and climate will be described. The present study area is located in the Western Dessert of Egypt as illustrated on figure 3.11 on page 19. As seen on the map, the area of interest consists of two oasis out of several. The cause of the location of the oasis is the groundwater in the area. The inhabited area in the dessert has increased during the last 50 years due to dessert reclamation, and it is still part of the national water resources policy to increase the area reclaimed from the dessert (Amer et al., 2005, p. 13). Construction starting in 1997, a megaproject called the Sheikh Zayed Canal will drain water from Lake Nasser at the Aswan dam and make it flow 850km through the dessert into the New Valley. This is expected to reclaim further at least 540.000 feddans1 of land (Ibrahim and Ibrahim, 2003, p. 252).
2.1
Geology
The groundwater underlying the Kharga Oasis and Dakhla Oasis is part of a larger aquifer, as Figure 2.1: Illustrative map illustrated by figure 2.2, known as the Nubian of the Nile Valley and Sandstone Aquifer System. The water in this New Valley, Egypt. Based aquifer stems from a period of a wetter climate, on (Waterbury and Whitand ”there is general agreement that the aquifer tington, 1998) is under a non-steady condition since the beginning of the development in Egyptian and Libyan oasis in 1960”(Gossel et al., 2004), meaning that the current discharge is larger than the recharge. However, the size of the aquifer (estimated 150.000km3 (Tahlawi et al., 2008)) means that the water can be extracted economically for many years to come. The current pumping rates are more than 3 million m3 /day in the Qattara Depression, 11
feddan = 4200m2
4
Figure 2.2: Boundary of the Nubian Sandstone Aquifer System (Grenmillion, 2010) 400.000 m3 /day in the Farafra-Baharyia Depression planning to increase to a total of 1 million m3 /day, and 3 million m3 /day in the Kharga-Dakhla Depression (Ibrahim and Ibrahim, 2003, p. 69). The Nubian Sandstone Aquifer is a confined aquifer with the main flow direction towards northeast. The aquifer is bounded to the north by the salinefreshwater borderline which runs beneath the Qattara Depression. To the south it is bounded by the Uweinat-Safsaf uplift and the Safsaf-Aswan uplift except between Gebel Kamil and the Bir Tafawi area where a south trending graben ensures a connection to Sudan (Thorweihe, 1990). To the east the aquifer is bounded by the the crystalline basement complex of the Nile as illustrated on figure 2.3. The aquifer is not bounded to the west, and therefore flow from the Kufra basin in Libya can occur. All of this means, that very little recharge of the aquifer, if any, occurs. The speed of the water is only 3cm/year through the Misha Graben from Sudan thus it would take 500 years to travel the distance of 1500km from the wet areas in Sudan (Ibrahim and Ibrahim, 2003, p. 45).The only interaction with the Nile occurs at Lake Nasser (Ebraheem et al., 2004).
5
Figure 2.3: Geological map of the Egyptian part of the Nubian Sandstone Aquifer (Ibrahim and Ibrahim, 2003, p. 43)
6
Figure 2.4: Geological transects of the Western Dessert of Egypt (Thorweihe, 1990)
7
Figure 2.5: Geological transects of the Western Dessert of Egypt (Thorweihe, 1990)
8
9
Name of Aquifer:
Type locality
Nubian Sandstone
Western Dessert Kharga Dakhla Bahariya Farafra East Oweinat
Depth of top aquifer(m)
Saturated thickness(m)
Depth to water table(m)
Hydraulic conductivity (m/day)
Porosity (%)
Salinity (ppm)
50–200 200 150–300 200–500 10–20
500–700 500–1000 1000–1500 1000–2000 100–300
0-30 0–20 0–20 Flow 20–30
2–4 6–7 5–10 2–5 10–20
20 20–25 20–40 20–30 20–30
< 1000 < 1000 < 1000 < 1000 < 1000
Table 2.1: Excerpt of Table of main aquifer systems of Egypt (Tahlawi et al., 2008)
The main characteristics of the aquifer in the study area are presented in table 2.1 on page 9. As seen from the table, the depth of the groundwater is increasing towards the north which makes it more expencive to recover the water. Moreover, the salt concentration is low over the whole area meaning that problems of salinization in agriculture are reduced.
2.2
Climate
As can be seen from figure 2.6 and figure 2.7 on page 11, there is a large variation in average temperature between summer and winther in the region. Since water so far has been provided free of charge, the amount of water used for irrigation has been adjusted to this fact. This is also confirmed by (Ali2 personal communication), stating that the evaporation pr. square meter is 16mm/day in July and 4mm/day in December. Other sources (Ahmed, 1999, p. 61) state the average evaporation rate to 18.4mm/day. The high average temperatures also mean that agriculture takes place all year round. As can also be seen from the figures, very little and very infrequent rain is falling in the area meaning that the all the agriculture is irrigated. This chapter has set the current situation in the study area into a geologic and climatic context. This serves as a background understanding for the results presented later in this report.
2 Dr.
Mohsen A. Gameh Ali, Assiyut University, Egypt
10
Figure 2.6: Temperature and precipitation in El-Kharga during the period of study. Data source (National Climatic Data Center, National Oceanic and Atmospheric Administration, United States Department of Commerce., 2011)
Figure 2.7: Temperature and precipitation in El-Dakhla during the period of study. Data source (National Climatic Data Center, National Oceanic and Atmospheric Administration, United States Department of Commerce., 2011) 11
Chapter 3
Method To estimate the temporal influence on the aquifer underlying Kharga Oasis and Dakhla Oasis in the Western Dessert of Egypt an analysis of data from the Gravity Recovery And Climate Experiment satellite remote sensing mission has been performed. This is chosen since obtaining in-situ groundwater measurements in a developing country like Egypt is very difficult. Moreover, data obtained from remote sensing sources are often cheaper than in-situ measurements, and it is possible through remote sensing data to obtain spatially continous measurements as opposed to point measurements obtained in-situ. Strassberg, Scanlon, and Rodell (Strassberg et al.) compared the amount of seasonal groundwater measured from GRACE, subtracted the soil moistrue from the NLDAS model, over the High Plains aquifer in USA, with in situ measurements of groundwater level. The High plains aquifer is an unconfined aquifer as opposed to the Nubian Sandstone Aquifer. The study found good agreement between calculated groundwater level and measured groundwater level. Rodell et al. (2007) examined groundwater changes in the Mississippi river basin, using GRACE and GLDAS, and compared to well measurements in the unconfined aquifers of the basin. The results show good agreement between satellite measurements and in-situ measurements for large basins ∼900.000km2 , but poorer results for basins∼500.000km2 . Tang et al. (2010) investigated the Klamath and Sacramento river basins in USA to examine the seasonal variations in terrestrial water storage (TWS) from GRACE compared to a combination of satellite data and surface data. The results show that GRACE substantially underestimate the seasonal cycle. As can be seen from this, by no means complete review of earlier GRACE studies, GRACE is a developing method in hydrological studies. Therefore this study also contributes to the knowledge about GRACE in relation to studies in arid regions and studies on small spatial scales. A total of 102 temporally unevenly spaced GRACE datasets with a temporal resolution of approximately 30 days were obtained from Center for Space Research (CSR) at the University of Texas, Austin covering the period from 2002 to 2010. Each dataset consist of a set of Stokes coefficients, representing the change in gravity field compared to the GGM02 gravity field model (Tapley et al., 2005), to degree and order 100. The fact that GRACE resolves to degree and order 100, means that the time-variable gravity changes are limited to struc12
Figure 3.1: Gaussian filter used to suppress errors in the data. tures of about 300km half wavelength (Schmidt et al., 2008). The used datasets are release 04 using a time averaged background model with a non-tidal atmosphere and ocean, which according to (Andersen1 , personal communication) is suitable for the present purpose. The obtained data are processed to remove the effect of the atmosphere (Flechtner, 2007). The temporal change in gravity was calculated from the Stokes coefficients as described in chapter 4. The data was smoothed using a Gaussian filter of halfwidth 0km in order to minimize the signal leakage and to suppress non-hydrological errors, and a spatial average were taken to improve the accuracy. The code of the implementation can be found in appendix A. The Gaussian filter was validated by calculating a large number of Stokes coefficients and through visual inspection noting that it converged towards the delimitation matrix. The calculated raster data were thus exported to GIS to map with data from other sources. The raster data were then validated quantitatively and qualitatively against raster data obtained from The Space Geodesy Research Group (GRGS)2 , as recommended by (Andersen, personal communication). GRGS has a different data processing procedure compared to this study (Horwath et al., 2011). The qualitative comparison of the calculated raster data and the raster data obtained from GRGS, show some resemblance in the area around Egypt as illustrated by figure 3.2 and figure 3.3. The large signal in the Mediterranean sea is located in the same spot om both images. Also comparing the global data on figure 3.3 and figure 3.4 shows some resemblance over e.g. India, South America and, the South Pole. Also looking at figure 3.4 alone, it is possible to see the 1 Senior
researcher at DTU space Ole B. Andersen
2 http://grgs.omp.obs-mip.fr/
13
Figure 3.2: Localized calculation of the geopotenital around Egypt for the period December 2010 based on data obtained from CSR continents clearly, indicating that the implementation is correct. However, the GRGS data are heavily influenced by striping, which is a satellite artefact. This is because the GRGS data are integrated over a ten day period, whereas the CSR data are integrated over 30 day period. To see if the images were actually correlated the correlation between the corresponding dataset for the period 2004 to 2010, where the data are continous was calculated, and the results shown on figure 3.5. The monthly data based on CSR were analyzed against the three corresponding dataset from GRGS. As can be seen on figure 3.5, the correlation between the datasets is very low. This can also be seen on figure 3.7, where an example of a scatterplot of data obtained from GRGS vs. data calculated from CSR is plotted. The reason for the low correlation is probably due to the short time period covered by the data from GRGS. As can be seen on figure 3.6, the striping effect is dominating, and thus the internal variance in the dataset is probably larger than any correlation. In this way an eventually small correlation is covered by a large internal variance. Nevertheless, the qualitative comparison indicates that the data calculated from CSR are valid. In order to delimit the hydrological region of interest, a map of the location of the oasis in the Western Dessert was digitized as illustrated on figure 3.8. This was done through scanning the map, and georeferencing by hand in GIS. A total of 17 ground control points was pointed out on the scanned map, and rectified in relation to vector data obtained from other sources. Next a second degree polynomial fit was performed to rectify the scanned image. The second degree polynomial fit was chosen since it yielded a better RMS, and because it by visual inspection yielded a better correspondance between the straight lines on the scanned map
14
Figure 3.3: Global raster representation of the geopotenial from day 348-357 (source GRGS),
Figure 3.4: Global raster data set for the period December 2010 calulated on the basis of data from CSR.
15
Figure 3.5: Statistical correlation between the corresponding dataset from GRGS and dataset calculated from CSR for the period 2004-2010.
Figure 3.6: Raster data from BGI in a period with high correlation with data calculated from CSR 16
Figure 3.7: Scatterplot of data obtained from GRGS vs. data calculated from CSR from day 244-253 in 2005, where the correlation is supposed to be slightly higher.
Figure 3.8: Digitized map with the relevant oasis included. (source: Thorweihe (1990)) 17
Figure 3.9: GRACE data from day 121 to 139 in 2002 calculated from CSR and in the vector data. The resulting digitalization has a RMS of 0.06 degrees corresponding to approximately 7 km, which is a high value(Booth et al., 2002). However, considering the low spatial resolution of GRACE and the original geometry of the map, it must be considered acceptable. Next the calculated raster data based on CSR was evaluated through visual comparison with the digitized map to determine the hydrological area of interest. An example of the maps are shown in figure 3.9 and 3.10. On figure 3.9 it seems like there is a signal coming from outside Egypt in the south-western part of the map. In figure 3.10 there are indications of a signal coming from the Nile and/or the Red Sea. Thus the smallest possible delimitation of the area of interest, considering the resolution of GRACE has been chosen as illustrated on figure 3.11 Lastly a spatial average over the defined area was taken as described in chapter 4, and the results plotted as a time-series. A non-linear least squares fit was made to a formula of the following kind: ∆σ region (t) = C + αt + Acos(ωt + φ)
(3.1)
Where, ∆σ region is the time dependent change in surface density, C is a constant surface density, α is the slope of an eventual linear trend, A is the amplitude of an eventual seasonal oscillation, ω is the angular velocity defined as 2π/T where T is the period of the oscillation and φ is a phase angle representing how much the oscillation is displaced according to a pure cosine-function.
18
Figure 3.10: GRACE data from day 213 to 243 in 2002 calculated from CSR
Figure 3.11: Hydrological area of interest
19
Chapter 4
Theory This chapter will explain the fundamentals with respect gravitational measurement of groundwater. There will in this chapter be an exposition of how groundwater can be measured as temporal changes in the earth gravitational field, and what mathematical filtering techniques are used in the subsequent empiri chapter. The section will from newtonian mechanics and derivations through potential theory, explain how the earth gravitational potential expanded as a sum of spherical harmonics, is related to the surface mass density of the earth. The assumption is that on timescales relevant to the GRACE mission, the only change in surface mass density is resulting from movement of water masses in the earth surface.
4.1
Technical description of the GRACE mission
The GRACE mission consists of two satellites chasing each other in low earth near polar orbit. The distance from the earth surface to the satellites will initially be ∼500km, which is low compared to other satellites used for earth observation. This means that the GRACE satellites will experience more air drag. When the first satellite experiences an increased gravity, it will accelerate due to Newtonian mechanics of circular motion. This will increase the distance (initially 220km) between the satellites. This increase in distance measured using a K-band radar (a horn on the one satellite to send the radar pulse, and a reflector on the other satellite to reflect the pulse back again). On the same time both satellites are equipped with GPS-receivers, to measure their exact position and accelerometers to distinguish the acceleration resulting from gravity from the acceleration resulting from air drag. After a little while the second satellite will measure the increase in gravity and speed up, thus decreasing the distance again. This distance change leads to a calculation of the geopotential as described below. Moreover both satellites are equipped with star-cameras and gyros to handle the internal allignment of the satellites.
20
Figure 4.1: GRACE satellites with main instruments. KBR-SST is the Kband Range Satellite-to-Satellite Tracking and GPS-SST is the GPS Satelliteto-satellite tracking (Schmidt et al., 2008)
Figure 4.2: The spatial and temporal resolution og GRACE. As can be seen there is a tradeoff between spatial and temporal resolution in the GRACE mission (Schmidt et al., 2008)
21
4.2
Introduction to potential theory
The GRACE satellite measures the temporal change in the gravity field of the earth. According to Newtonian mechanics, the gravitational force can be described in vector formulation: K = −G m1l2m2 ll Where K is the gravitational force, G is the gravitational constant, m1 is the mass of the object subjected to attraction, m2 is the the mass of the attracting object and l is the distance between the center of mass of the two objects. Setting the mass of the attracted object to unity transforms the gravitational force into the gravitational attraction(b): l b = −G m l2 l
The earth gravitational field is easier described through the scalar quantity the potential (V ), rather than the vector quantity the acceleration. Since ∇×K=0 the gravitational force can be described as the gradient of the potential: b = ∇V For the earth approximated as a point mass the potential will be(Torge, 2001, p. 45–47): V = GM (4.1) l ZZZ ρ V (r) = G l dv earth
By applying the law of cosines to the OP P ′ triangle as illustrated in figure 4.3
Figure 4.3: The gravitational coordinate system with the attracting mass located in P ′ and the attracted mass located in P , and the origo of the coordinate system in the earth center of mass.
22
we obtain: 1 l
If
1 l
2
′2
′
= r + r − 2rr cos(ψ)
− 1
2
=
1 r
1+
′ 2 r r
−
′ 2 rr cos(ψ)
− 1
2
is expanded into a power series, the following result is obtained: 1 l
=
1 r
∞ l X ′ r r
Pl (cos(ψ))
(4.2)
l=0
Where Pl (cos(ψ)) represent the l’th degree Legendre polynomial in cos(ψ). A decomposition of Pl (see e.g. Heiskanen and Moritz 1967, p. 33) results in: Pl (cos(ψ)) = Pl (cos(θ))Pl (cos(θ′ ) +2
l X (1−m)!
(1+m)! (Plm (cos(θ)cos(mλ)Plm (cos(θ
′
)cos(mλ′ )
m=1
+ Plm (cos(θ))sin(mλ)Plm (cos(θ′ )sin(mλ′ )
(4.3)
Where θ is the colatitude of the earth and λ is the longitude of the earth. Inserting equation 4.2 and 4.3 into equation 4.1 yields the following expression for the gravitational potential(Torge, 2001, p. 67–70): ! l ∞ X X a l GM 1+ (Clm cos(mλ) + Slm sin(mλ)) Plm (cos(θ)) (4.4) V = r r L=1M =0
Where a is the average radius of the earth.
4.3
The disturbing potential
When measuring the gravity field of the earth from a satellite, the earth can no longer be treated as a perfect sphere. In reality the earth density is inhomogeneous and position dependent, and so is the gravity field. The simplest approximation of the earths gravity field is an ellipsoid(U ) called the normal gravity potential(Torge, 2001, p. 103). The normal gravity potential is furthermore required to coincide as close as possible with measurements of the actual gravity field(W ). The difference between the normal actual gravity potential and the normal gravity potential is called the disturbing potential(T ). The disturbing potential in a point can thus be expressed as(Torge, 2001, p. 214): Tp = Wp − U p
(4.5)
Since the disturbing potential is a harmonic function and thus has to satisfy the Laplace equation, it can be expressed in the form of spherical coordinates: Tp =
GM r
l ∞ X X
L=1M =0
a l r
(Clm cos(mλ) + Slm sin(mλ)) Plm (cos(θ))
(4.6)
This is the form of the GRACE data, since cf. section 3, the GRACE data are gravity potential deviations from the GGM02 model. This model is constructed 23
with a constant potential (U ) and the same geometric shape as the physical earth called the telluriod(Torge, 2001, p. 112). The geoid used in the modelling also has to satisfy Laplace equation. A Taylor development of the normal potential in the telluroid point Q results in(Torge, 2001, p. 257): Up = UQ + ∂U ∂n ζp + . . .
Where n is the normal to U = UQ and the normal gravity defined as: γQ = − ∂U ∂n Results in the height anomaly(ζ):
ζp =
Tp −(Wp −UQ ) γQ
U Q = WP This means, that the geoid height anomaly can be expressed as: ζp =
Tp γp
(4.7)
Thus the geoid height anomaly can be expressed in spherical coordinates as: ζ(r, θ, λ) =
GM rγ
l ∞ X X
L=1M =0
a l r
(Clm cos(mλ) + Slm sin(mλ)) Plm (cos(θ))
(4.8)
Close to the earth surface it is possible to approximate the normal gravity (because the difference in axes is very small and small compared to the scale of GRACE) as a sphere such that: γ = GM r2 r=a
(4.9) (4.10)
This is also supported as the first term in a Taylor expansion. Inserting this into equation 4.8 yields the following expression for the geoid height anomaly. ζ(r, θ, λ) = r
l ∞ X X
(Clm cos(mλ) + Slm sin(mλ)) Plm (cos(θ))
(4.11)
L=1M =0
4.4
Relationship between the stokes coefficients and the earth surface density
This section is based on Chao and Gross 1987. The aim of this section is to explain how the geopotential can be used to probe the earth crust. This is done, by relating the geopotential as expressed in equation 4.11 to the earth surface density defined as the mass pr. area. By integrating over the selected area it is possible to determine the change in mass over time and thus also in volume. Expanding the geopotential as in equation 4.4 on page 23 as a power series and a multipole expansion yields the following relation between the Stokes coefficients and the multipole moment of the density distrubution (Γ): Clm + iSlm = 24
Klm Γ M al lm
(4.12)
Where Clm and Slm are the Stokes coefficients, i is the imaginary number, M is the mass of the earth, a is a scale factor, l is the Legendre degree and m is the Legendre order. The multipole moment of the density distribution is defined as: Z Γlm = ρ(r)rl Ylm (Ω)dV (4.13) Where ρ(r) is the density distribution function, r is the average radius of the earth, l is the Legendre degree, m is the Legendre order and Y is the surface spherical harmonic functions as a function of the angles abbreviated by Ω. The surface spherical harmonic functions are defined as: Ylm (Ω) = (−1)m
(2l+1)(l−m)! 4π(l+m)!
1
2
Plm (cos(θ))eimφ
(4.14)
Where Plm is the associated Legendre function. Inserting equation 4.13 and 4.14 into equation 4.12 yields the following relationship between the Stokes coefficients and the surface density of the earth (Wahr et al., 1998): Z l+2 cos(mφ) ∆Clm sin(θ)dθdφdr = 4πaρave3 (2l+1) ∆ρ(r, θ, φ)Pelm (cos(θ)) ar ∆Slm sin(mφ) (4.15) The following section is based on Wahr et al. 1998. Since the layers in the surface of the earth, where the mass is changing within a timescale relevant to the GRACE mission are very thin compared to the radius of the earth, a surface density can be defined as the radial integral through this layer: Z ∆ρ(r, θ, φ)dr (4.16) σ(θ, φ) = thin layer
l+2 ≈ The above mentioned approximation simplifies equation 4.16 because ar 1: Z cos(mφ) ∆Clm 3 e sin(θ)dθdφdr = 4πaρave (2l+1) ∆ρ(r, θ, φ)Plm (cos(θ)) sin(mφ) ∆Slm surface mass (4.17) The temporal redistribution of surface mass causes loading and deformation of the underlying solid earth. This causes an increased contribution to the Stokes coefficients as a factor kl called the load Love number of degree l. Z ∆Clm 3kl elm (cos(θ)) cos(mφ) sin(θ)dθdφdr ∆ρ(r, θ, φ) P = 4πaρave (2l+1) sin(mφ) ∆Slm solid earth (4.18) The load Love number is calculated based on a linear interpolation to the data presented in Wahr et al. 1998, which will result in errors of less than 0.05% (Wahr et al., 1998). The total change in geoid must then be the sum of the contribution from surface mass redistribution and elastic deformation of the solid earth: ∆Clm ∆Clm ∆Clm = + (4.19) ∆Slm ∆Slm surface mass ∆Slm solid earth
25
Since the temporal change in surface density has to obey the Laplace equation it can be decomposed in the following way: ∆σ(θ, φ) = aρw
l ∞ X X
L=1M =0
elm cos(mλ) + Selm sin(mλ) Plm (cos(θ)) C
(4.20)
elm and Selm are the Stokes coefficients Where a is the average earth radius, C of both the surface mass redistibution and the elastic deformation of the solid earth. From equation 4.17 and 4.18 the following relationship between the e S) e and the Stokes coefficients resulting from corrected Stokes coefficients (C, the GRACE measurements: ( ) elm ∆Clm ∆C ρave 2l+1 (4.21) = 3ρw 1+kl ∆Slm ∆Selm
Inserting equation 4.21 into equation 4.20 yields the following relationship between the earth surface density and the Stokes coefficients obtained from the GRACE mission. ∆σ(θ, φ) =
aρave 3
l ∞ X X
l=0 M =0
4.5
2l+1 (∆Clm cos(mφ) + ∆Slm sin(mφ)) Pelm (cos(θ)) 1+k l
(4.22)
Handling of GRACE error
The data obtained from the GRACE mission will contain errors resulting from the satelite instrumentation and from the truncation of the infinite series at a certain Legendre degree. As seen from figure 4.4, the satellite errors are getting larger proportionally to the degrees. The aim of the filtering will therefore be to reduce the amplitude of the higher degrees, thus reducng the satellite errors at the loss of spatial resolution. The errors can be reduced by spatially averaging the data leading to a better representation of the gravitational potential at the expence of spatial resolution.This section is based on Swenson and Wahr 2002. Setting up a function delimiting the region of interest according to the following: 0 Outside the basin ϑ(θ, φ) = (4.23) 1 Inside the basin The spatial average over the region is then calculated as: Z 1 ∆σ(θ, φ)ϑ(θ, φ)dΩ ∆σ region = Ωregion
(4.24)
Expanding the spatial average as a sum of spherical harmonics yields: ∆σ region =
aρave 3Ωregion
∞ X l X
2l+1 1+kl
(ϑclm ∆Clm + ϑslm ∆Slm )
(4.25)
l=0 M =0
Where ϑclm and ϑslm are the spherical harmonic coefficients representing ϑ(θ, φ): ϑ(θ, φ) =
1 4π
∞ X l X
l=0 M =0
Pelm (cos(θ)) (ϑclm cos(mφ) + ϑslm sin(mφ)) 26
(4.26)
Figure 4.4: Satellite measurement errors as a function of degree.
ϑclm ϑslm
=
Z
ϑ(θ, φ)Pelm (cos(θ))
cos(mφ) sin(mφ)
dΩ
(4.27)
If instead of using the formula presented in equation 4.23, using a spatial average(W ) to delimit the area: Z f region = 1 ∆σ(θ, φ)W (θ, φ)dΩ (4.28) ∆σ Ωregion Expanding the spatial average(W ) as a sum of spherical harmonics yields: W (θ, ω) =
1 4π
l ∞ X X
l=0 M =0
c s Pelm (cos(θ)) (Wlm cos(mφ) + Wlm sin(mφ))
This means that the approximate spatial average can be expressed as: X Kl c s f region = ∆σ Ωregion (Wlm ∆Clm + Wlm ∆Slm )
(4.29)
(4.30)
l,m
Kl =
aρE (2l+1) 3 1+kl
(4.31)
The definition of a Gaussian filter between two points represented by the angle between the two points (γ)(Jekeli, 1981): W (γ) = b=
b e−b(1−cos(γ)) 2π 1−e−2b ln(2) ! 1−cos r 1 /a 2
27
(4.32) (4.33)
Where r 1 /a is the half width of the Gaussian smoothing function. The Gaussian 2
smoothing function changes smoothly from a value of 1 inside the boundary to a value of 0 outside the boundary over the distance of approximately r 1 . The weighing coefficients are calculated as follows: c c θlm Wlm = 2πW l s s θlm Wlm
2
(4.34)
In this section the mathematical foundations of the data processing has been derived. These are the results implemented in Matlab as described in 3, and the results are presented in the following section.
28
Chapter 5
Results and conclusion The results of the project is illustrated in figure 5.1, along with the least squares fit to the function. The fit parameters are listed in table 5.1. As can be seen from both the figure and the table, GRACE measures a strong annual signal. This is also seen when calculating the oscillatory frequency into a number of days: 2π T = = 365.3300days 0.0172days−1 When seing that the amount of water is smallest in the late summer, this could indicate that it is a result of irrigation. If this is the case, it means that a large amount of water is pumped during the summer, and when the temperature drops, the irrigation is turned down, meaning that the groundwater has time to run into the area. The amplitude of 41.3142 Kg m2 is large compared to the earlier stated 3 million cubic meters of pumping. This could however be because the stated amount of pumping is erroneous, like many official numbers in a developing country like Egypt. The linear trend is very small, and thus it is more likely a result of satellite measurement errors than an actual hydrological effect.
Variable: C α A ω φ
Value: -1.1318 Kg m2 0.0010 m2Kg ·day 41.3142 Kg m2 0.0172 days−1 1.3891
Table 5.1: Fit parameters for the least squares fit to data. 29
Figure 5.1: Temporal development in terrestrial water content in the years 20022010.
30
Bibliography Adriansen, H. K. (2009). Land reclamation in Egypt: A study of life in the new lands. Geoforum 40, 664–674. Ahmed, A. A.-S. (1999). Geoelectrical and Hydrogeological Studies on the Groundwater Resources in Bulaq – Baris area, Kharga Oasis, Egypt. Master’s thesis, Assiut University. Amer, M. H., M. B. Attia, H. Fahmy, and M. Tawfik (2005). Water Policy Issues of Egypt. Technical report, International Commission on Irrigation and Drainage. Booth, B., J. Shaner, S. Crosier, P. Sanchez, and A. MacDonald (2002). Editing in ArcMap. ESRI. http://www.zmuc.dk/public/GIS-intro/Litteratur/ ESRIdigitalBooks/Editing_in_ArcMap.pdf. Chao, B. F. and R. S. Gross (1987). Changes in Earth’s rotation and low-degree gravitational field induced by earthquakes. Geophysical journal of the royal astronomical society 91, 569–596. Ebraheem, A. M., S. Riad, P. Wycisk, and A. M. Sefelnasr (2004). A local-scale groundwater flow model for groundwater resources management in Dakhla oasis, SW Egypt. Hydrogeology Journal 12, 714–722. Flechtner, F. (2007). AOD1B Product Description Document for Product Releases 01 to 04. Technical report, GeoForschungszentrum Potsdam, Department 1: Geodesy and Remote Sensing. http://op.gfz-potsdam.de/grace/ results/grav/AOD1B_PDD_20070413.pdf. Gossel, W., A. M. Ebraheem, and P. Wycisk (2004). A very large scale GISbased groundwater flow model for the nubian sandstone aquifer in Eastern Sahara (Egypt, northern Sudan and eastern Libya). Hydrogeology Journal 12, 698–713. Grenmillion, P. (2010). New Light Shed on the Nubian Aquifer. Water & Environment News (26), 2–5. Heiskanen, W. A. and H. Moritz (1967). Physical Geodesy. W. H. Freeman and Company. Horwath, M., J.-M. Lemoine, R. Biancale, and S. Bourgogne (2011). Improved GRACE Science results after adjustment of geometric biases in the Level-1B K-band ranging data. Journal of Geodesy 85, 23–38. 31
Ibrahim, F. N. and B. Ibrahim (2003). I.B.Tauris & Co. Ltd.
Egypt: An Economic Geography.
Jekeli, C. (1981). Alternative methods to smooth the Earth’s gravity field. Technical report, Department of Civil and Environmental Engineering and Geodetic Sciences, The Ohio State University, Columbus. http://www.geology. osu.edu/~jekeli.1/OSUReports/reports/report_327.pdf. Meyer, G. (2004). Die Arabische Welt Im Spiegel der Kulturgeographie, Chapter Neues Bew¨asserungsland in der ¨agyptischen W¨ uste: von Kleinbauern und Staatsbetrieben zu freizeitfarmern und Großinvestoren. Zentrum f¨ ur Forschung zur Arabische Welt. Mitchell, T. (1995). The Object of Development: America’s Egypt, Chapter 7, pp. 125–151. Routledge. National Climatic Data Center, National Oceanic and Atmospheric Administration, United States Department of Commerce. (2011). http://www.ncdc. noaa.gov. Popp, H. (2004). Die Arabische Welt Im Spiegel der Kulturgeographie, Chapter Einf¨ uhrung. Zentrum f¨ ur Forschung zur Arabische Welt. Rodell, M., J. Chen, H. Kato, J. S. Famiglietti, J. Nigro, and C. R. Wilson (2007). Estimating groundwater storage changes in the Mississippi River basin (USA) using GRACE. Hydrogeology Journal 15, 159–166. Schmidt, R., F. Flechtner, U. Meyer, K.-H. Neumayer, C. Dahle, R. K¨ onig, and J. Kusche (2008). Hydrological Signals Observed by the GRACE satellites. Survey of Geophysics 29, 319–334. Strassberg, G., B. R. Scanlon, and M. Rodell. Comparison of seasonal terrestrial water storage variations from GRACE with groundwater-level measurements from the H]igh Plains aquifer, journal = Geophysical Research Letters, year = 2007, volume = 34, pages = LI4402,. Swenson, S. and J. Wahr (2002). Methods for inferring regional surface-mass anomalies from Gravity Recovery and Climate Experiment(GRACE) measurements of time-variable gravity. Journal of geophysical research 107 (B9), 2193. Tahlawi, M. R. E., A. A. Farrag, and S. S. Ahmed (2008). Groundwater of Egypt: ”an environmental overview”. Environmental Geology 55, 639–652. Tang, Q., H. Gao, P. Yeh, T. Oki, F. Su, and D. P. Lettenmaier (2010, February). Dynamics of Terrestrial Water Storage Change from Satellite and Surface Observations and Modeling. Hournal of Hydrometeorology 11, 156–170. Tapley, B., J. Ries, S. Bettadpur, D. Chambers, M. Cheng, F. Condi, B. Gunter, Z. Kang, P. Nagel, R. Pastor, T. Pekker, S. Poole, and F. Wang (2005). GGM02 – An improved Earth gravity field model from GRACE. Journal of Geodesy 79, 467–478. Thorweihe, U. (1990). The Geology of Egypt, Chapter Nubian Aquifer system, pp. 601–611. A.A.Balkema/Rotterdam/Brookfield. 32
Torge, W. (2001). Geodesy (3rd edition ed.). Walter de Gruyter. Wahr, J., M. Molenaar, and F. Bryan (1998). Time variabillity of the Earth’s gravity field: Hydrological and oceanic effects and their possible detection using GRACE. Journal of geophysical research 103 (B12), 30205–30229. Waterbury, J. and D. Whittington (1998). Playing chicken on the Nile? The implications of microdam development in the Ethiopian highlands and Egypt’s New Valley Project. Natural Resources Forum 22 (3), 155–163. World Commission on Environment and Development (1987). Our Common Future. Oxford University Press.
33
Appendix A
Implementation in MatLab Listing A.1: Code for calculating the geopotential, the Gaussian filter and the spatial average 1
% Program t o c a l c u l a t e t h e t e m p o r a l change i n t h e g e o p o t e n t i a l from t h e % GRACE s a t e l l i t e s .
3 5 7
% Beregner Load Love Numbers : f or i = 1 : 2 0 1 Kl3 = interp1 ( l , Kl , Kl2 ) ; %L i n e a r i n t e r p o l a t i o n bet w een t h e g i v e n d a t a p o i n t s end
9 11
13 15
17
19 21 23
% The Gaussian f i l t e r i s r e s e t VarTheta = zeros ( round ( ( EndLat−S t a r t L a t ) / D e l t a ) +1,round ( ( EndLong−St ar tL ong ) / D e l t a ) +1) ; f or l =1:round ( ( EndLat−S t a r t L a t ) / D e l t a )+1 l a t=round ( ( EndLat−S t a r t L a t ) / D e l t a )+2− l ; P1 = 90−( S t a r t L a t +( l −1)∗ D e l t a ) ; %C o n v e r t i n g from l a t i t u d e t o c o l a t i t u d e % The g e o g r a p h i c a l c o o r d i n a t e o f t h e l a t i t u d e i s converted to radians : pos1 = 2∗ pi ∗ ( P1 / 3 6 0 ) ; f or m=1:round ( ( EndLong−St ar tL on g ) / D e l t a )+1 l o n g = m; P2 = Sta rt Lo ng +(m−1)∗ D e l t a ; % The g e o g r a p h i c a l l o n g i t u d e i s t r a n s f o r m e d t o radians : pos2 = 2∗ pi ∗ ( P2 / 3 6 0 ) ;
25
%
Calculation of the geopotential 34
Temp ( : , 1 ) = DataReal ( : , 1 , nr ) ; % Degree nr i s s t o r e d Temp ( : , 2 ) = DataReal ( : , 2 , nr ) ; % Order nr i s s t o r e d Temp ( : , 3 ) = ( DataReal ( : , 3 , nr ) . ∗ cos ( DataReal ( : , 2 , nr ) . ∗ pos2 ) ) +(DataReal ( : , 4 , nr ) . ∗ sin ( DataReal ( : , 2 , nr ) ∗ pos2 ) ) ; Temp ( : , 3 ) = Temp ( : , 3 ) . ∗ ( ( 2 ∗ DataReal ( : , 1 , nr ) +1) ./(1+ KlReal ( ( DataReal ( : , 1 , nr ) +1) , 2 ) ) ) ; Temp = s o r t r o w s (Temp) ; % R e o r d e r i n g a c c o r d i n g t o d e g r e e and o r d e r Temp ( : , 3 ) = Temp ( : , 3 ) . ∗ Plm2 ( : , 3 + l −1) ; % M u l t i p l i c a t i o n with f u l l y normalized legendre polynomial Sigma ( l a t , l o n g ) = sum(Temp ( : , 3 ) ) ; % Adding up t h e r e s u l t Sigma ( l a t , l o n g ) = Re∗RhoE/3 ∗ Sigma ( l a t , l o n g ) ;
27
29
31
33
35
% C a l c u l a t i o n o f t h e Gaussian f i l t e r ThetaTemp ( : , 1 ) = DataReal ( : , 1 , nr ) ; ThetaTemp ( : , 2 ) = DataReal ( : , 2 , nr ) ; ThetaTemp ( : , 3 ) = Plm2 ( : , 3 + l −1) . ∗ (WlmC( : , 3 ) . ∗ cos ( Plm2 ( : , 2 ) ∗ pos2 ) . . . +WlmS ( : , 3 ) . ∗ sin ( Plm2 ( : , 2 ) ∗ pos2 ) ); VarTheta ( l a t , l o n g ) = sum( ThetaTemp ( : , 3 ) ) ; VarTheta ( l a t , l o n g )=VarTheta ( l a t , l o n g ) + Plm2 ( i ,3+ l −1) ∗ (WlmC( i , 3 ) ∗ cos ( Plm2 ( i , 2 ) ∗ pos2 ) . . . +WlmS( i , 3 ) ∗ sin ( Plm2 ( i , 2 ) ∗ pos2 )); VarTheta ( l a t , l o n g )=VarTheta ( l a t , l o n g ) ∗ ( 1 / ( 4 ∗ pi ) ) ;
37 39
41
43
45
% Calculation of the f i l t e r e d geopotential AppSigmaTemp = s o r t r o w s ( DataReal ( : , : , nr ) ) ; AppSigma (R, nr ) = sum( ( ( 2 ∗ DataReal ( : , 1 , nr ) +1) ./(1+ KlReal ( ( DataReal ( : , 1 , nr ) +1) , 2 ) ) ) . . . . ∗ (WlmC( : , 3 ) . ∗ AppSigmaTemp ( : , 3 )+ WlmS ( : , 3 ) . ∗ AppSigmaTemp ( : , 4 ) ) ); AppSigma (R, nr ) = AppSigma (R, nr ) ∗ ( Re∗RhoE) / ( 3 ∗ Angular Area ) ;
47
49
end
51
end Listing A.2: Code for calculating the correlation of the two datasets. %−−−−−−−−−−−−−−−Program t o c a l c u l a t e t h e c o r r e l a t i o n bet w een two d a t a s e t −−−− 2
35
4 6
% F i r s t t h e mean v a l u e s a r e s u b t r a c t e d t o mode t h e c o o r d i n a t e system t o t h e % center of the dataset . GRACE = reshape (GRACEGlob ( : , : , nr ) , 1 8 0 ∗ 3 6 0 , 1 ) ; f or j = 1 : 3 BGICca = reshape ( BGIcc ( : , : , 4 0 + j −19) , 1 8 0 ∗ 3 6 0 , 1 ) ;
8 10 12
MiddelGRACE = mean(GRACE) ; MiddelBGI = mean( BGICca ) ; GRACE = GRACE − repmat (MiddelGRACE , 1 8 0 ∗ 3 6 0 , 1 ) ; BGICca = BGICca − repmat ( MiddelBGI , 1 8 0 ∗ 3 6 0 , 1 ) ;
14
16 18
% GRACE and BGICCa a r e merged and t h e v a r i a n c e −c o v a r i a n c e matrix i s % calculated XCov = cov ( [GRACE BGICca ] ) ; Rho ( nr−19+ j ) = XCov ( 1 , 2 ) / ( sqrt (XCov ( 1 , 1 ) ) ∗ sqrt (XCov ( 2 , 2 ) ) ); end Listing A.3: Code for calculating the coefficients used in the Gaussian filter
1
%Program t o c a l c u l a t e t h e c o e f f i c i e n t s used f o r t h e Gaussian f i l t e r i n t h e %d a t a i n t e r p r e t a t i o n o f GRACE.
3 5 7
f or n = 1 : 5 1 5 1 % Number o f S t o k e s c o e f f i c i e n t s f or l =1:round ( ( EndLat−S t a r t L a t ) / D e l t a )+1 l a t=round ( ( EndLat−S t a r t L a t ) / D e l t a )+2− l ; % L a t i t u d e i s c o n v e r t e d t o c o l a t i t u d e and t o r a d i a n s P1 = 90−( S t a r t L a t +( l −1)∗ D e l t a ) ; pos1 = 2∗ pi ∗ ( P1 / 3 6 0 ) ; f or m=1:round ( ( EndLong−St ar tL on g ) / D e l t a )+1 l o n g = m; P2 = Sta rt Lo ng +(m−1)∗ D e l t a ; % The l o n g i t u d e i s c o n v e r t e d t o r a d i a n s pos2 = 2∗ pi ∗ ( P2 / 3 6 0 ) ;
9 11 13 15 17
%I n t e g r a l o f t h e d a t a o v e r s o l i d a n g l e VarThetaC ( n , 3 ) =VarThetaC ( n , 3 )+A r e a D e l i m i t ( l a t ,m ) ∗ sin ( pos1 ) ∗ ( ( ( D e l t a ∗ 1 1 1 0 0 0 ) ˆ 2 ) / ( ( Re ) ˆ 2 ) ) ∗Plm2 ( n,3+ l −1)∗ cos ( Plm2 ( n , 2 ) ∗ pos2 ) ; VarThetaS ( n , 3 ) =VarThetaS ( n , 3 )+A r e a D e l i m i t ( l a t ,m ) ∗ sin ( pos1 ) ∗ ( ( ( D e l t a ∗ 1 1 1 0 0 0 ) ˆ 2 ) / ( ( Re ) ˆ 2 ) ) ∗Plm2 ( n,3+ l −1)∗ sin ( Plm2 ( n , 2 ) ∗ pos2 ) ;
19
end
21
end 36
% Degree and o r d e r numbers a r e s t o r e d VarThetaC ( n , 2 ) = Plm2 ( n , 2 ) ; VarThetaC ( n , 1 ) = Plm2 ( n , 1 ) ; VarThetaS ( n , 2 ) = Plm2 ( n , 2 ) ; VarThetaS ( n , 1 ) = Plm2 ( n , 1 ) ; WlmC( n , 2 ) = Plm2 ( n , 2 ) ; WlmC( n , 1 ) = Plm2 ( n , 1 ) ; WlmS( n , 2 ) = Plm2 ( n , 2 ) ; WlmS( n , 1 ) = Plm2 ( n , 1 ) ; % The f i r s t S t o k e s c o e f f i c i e n t i s c a l c u l a t e d i f Plm2 ( n , 1 )==0 count =1; Wl( n , 2 ) = 1 / ( 2 ∗ pi ) ; Wl( n , 1 )=Plm2 ( n , 1 ) ; end % The second S t o k e s c o e f f i c i e n t i s c a l c u l a t e d i f Plm2 ( n , 1 )==1 Wl( n , 2 ) = 1 / ( 2 ∗ pi ) ∗(((1+ exp(−2∗b ) ) /(1−exp(−2∗b ) ) ) − 1/b ) ; Wl( n , 1 )=Plm2 ( n , 1 ) ; end % S t o k e s c o e f f i c i e n t s o f h i g h e r d e g r e e and o r d e r a r e calculated i f Plm2 ( n , 1 ) >1 i f Plm2 ( n , 1 ) ˜=Plm2 ( n −1 ,1) Wl( n , 2 ) = −((2∗Plm2 ( n , 1 ) +1)/b ) ∗Wl( n −1 ,2)+Wl( count , 2 ) ; count=n−1; else Wl( n , 2 )=Wl( n −1 ,2) ; end Wl( n , 1 )=Plm2 ( n , 1 ) ; end WlmC( n , 3 ) = 2∗ pi ∗Wl( n , 2 ) ∗VarThetaC ( n , 3 ) ; WlmS( n , 3 ) = 2∗ pi ∗Wl( n , 2 ) ∗ VarThetaS ( n , 3 ) ;
23 25 27 29 31 33 35 37 39
41 43
45
47 49 51 53 55
end Listing A.4: Code for calculating the fully normalized associated spherical harmonics
1
3 5 7
%Program t o c a l c u l a t e t h e f u l l y n o r m a l i z e d a s s o c i a t e d s p h e r i c a l harmonics Plm2=zeros ( 5 1 5 1 , round ( ( EndLat−S t a r t L a t ) / D e l t a ) ) ; % Resetting the matrix m=2; f or i =1:round ( ( EndLat−S t a r t L a t ) / D e l t a )+1 P1 = 90−( S t a r t L a t +( i −1)∗ D e l t a ) ; pos1 = 2∗ pi ∗ ( P1 / 3 6 0 ) ; count =7; %The f i r s t s e v e n s p h e r i c a l harmonics a r e c a l c u l a t e d 37
9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39
41 43 45 47
49
manyally Plm2 (1 ,3+ i −1)=1; i f i ==1; Plm2 ( 1 , 1 ) =0; Plm2 ( 1 , 2 ) =0; end Plm2 (2 ,3+ i −1)=sqrt ( 3 ) ∗ cos ( pos1 ) ; i f i ==1 Plm2 ( 2 , 1 ) =1; Plm2 ( 2 , 2 ) =0; end Plm2 (3 ,3+ i −1)=sqrt ( 3 ) ∗ sin ( pos1 ) ; i f i ==1 Plm2 ( 3 , 1 ) =1; Plm2 ( 3 , 2 ) =1; end Plm2 (4 ,3+ i −1) =0.5∗ sqrt ( 5 ) ∗ ( 3 ∗ ( cos ( pos1 ) ) ˆ2−1) ; i f i == 1 Plm2 ( 4 , 1 ) =2; Plm2 ( 4 , 2 ) =0; end Plm2 (5 ,3+ i −1)=sqrt ( 1 5 ) ∗ sin ( pos1 ) ∗ cos ( pos1 ) ; i f i ==1 Plm2 ( 5 , 1 ) =2; Plm2 ( 5 , 2 ) =1; end Plm2 (6 ,3+ i −1) =0.5∗ sqrt ( 1 5 ) ∗ ( sin ( pos1 ) ) ˆ 2 ; i f i ==1 Plm2 ( 6 , 1 ) =2; Plm2 ( 6 , 2 ) =2; end f or l = 3 : 1 0 0 Plm2 ( count ,3+ i −1)=−(sqrt ( 2 ∗ l +1)/ l ) ∗ ( ( l −1) / ( sqrt ( 2 ∗ l −3) ) ) ∗Plm2 ( count −2∗ l +1,3+ i −1) . . . +(sqrt ( 2 ∗ l +1)/ l ) ∗ sqrt ( 2 ∗ l −1)∗ cos ( pos1 ) ∗Plm2 ( count−l ,3+ i −1) ; i f i == 1 Plm2 ( count , 1 )=l ; Plm2 ( count , 2 ) =0; end count=count +1; f or m=2: l −1 Plm2 ( count ,3+ i −1) =(((2∗ l +1) ∗ ( 2 ∗ l −1) ) / ( ( l +(m −1) ) ∗ ( l −(m−1) ) ) ) ˆ ( 0 . 5 ) . . . ∗ cos ( pos1 ) ∗Plm2 ( count−l ,3+ i −1) . . . −(((2∗ l +1) ∗ ( l +(m−1)−1) ∗ ( l −( m−1)−1) ) / ( ( 2 ∗ l −3) ∗ ( l +(m −1) ) ∗ ( l −(m−1) ) ) ) ˆ ( 0 . 5 ) 38
... ∗Plm2 ( count −2∗ l +1,3+ i −1) ; i f i == 1 Plm2 ( count , 1 )=l ; Plm2 ( count , 2 )=m−1; end count=count +1;
51 53 55
end Plm2 ( count ,3+ i −1)=sqrt ( 2 ∗ l +1)∗ cos ( pos1 ) ∗Plm2 ( count−l ,3+ i −1) ; i f i == 1 Plm2 ( count , 1 )=l ; Plm2 ( count , 2 )=l −1; end count=count +1; Plm2 ( count ,3+ i −1) =sqrt ( ( 2 ∗ l +1) / ( 2 ∗ l ) ) ∗ sqrt (1 −( cos ( pos1 ) ) ˆ 2 ) ∗Plm2 ( count−l −1,3+ i −1) ; i f i == 1 Plm2 ( count , 1 )=l ; Plm2 ( count , 2 )=l ; end count=count +1;
57
59 61 63
65 67 69
end end
39
