Nature & Wildlife

Wideband Autocorrelation Radiometry for Lake Icepack Thickness Measurement With Dry Snow Cover

Description
Wideband autocorrelation radiometry (WiBAR) is a recently developed microwave radiometric technique to measure the lake icepack or snowpack thickness. This technique offers a direct method to remotely measure the microwave propagation time difference
Published
of 5
3
Published
All materials on our website are shared by users. If you have any questions about copyright issues, please report us to resolve them. We are always happy to assist you.
Similar Documents
Share
Transcript
  This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. IEEE GEOSCIENCE AND REMOTE SENSING LETTERS 1 Wideband Autocorrelation Radiometry forLake Icepack Thickness MeasurementWith Dry Snow Cover Seyedmohammad Mousavi ,  Student Member, IEEE  , Roger De Roo ,  Member, IEEE  ,Kamal Sarabandi ,  Fellow, IEEE  , Anthony W. England ,  Fellow, IEEE   Abstract —Wideband autocorrelation radiometry (WiBAR) is arecently developed microwave radiometric technique to measurethe lake icepack or snowpack thickness. This technique offers adirect method to remotely measure the microwave propagationtime difference of multipath microwave emission from low-losslayered surfaces such as dry snowpack and freshwater lakeicepack. The microwave propagation time difference through thepack yields a measure of its vertical extent. However, the lakeicepack thickness measurement can be affected by the presenceof a dry snowpack, which introduces another multipath interfer-ence. We present a simple geophysical forward model consideringthis effect and derive the WiBAR system requirements needed tocorrectly measure the icepack thickness. An X-band instrumentfabricated from commercial-off-the-shelf (COTS) components isused to measure the thickness of fresh water lake ice at theUniversity of Michigan Biological Station. The WiBAR was ableto directly measure the icepack of about 36 cm with a snowpackof about 4 cm on top at incidence angle of 69.4 ◦ with an accuracyof 2 cm.  Index Terms —Ice, microwave radiometry, multilayered media,multipath interference, remote sensing, snow. I. I NTRODUCTION M ONITORING the ice thickness is crucial in freshwaterresource management [1], [2]. Ice thickness is alsoimportant in analyzing the pressure exerted to off-shore struc-tures such as wind farms [3]. In addition, estimation of lakeice thickness is essential for the safety of ice fishing and iceskating activities. The current method of ice thickness mea-surement is by drilling holes through the ice, which is not onlydifficult but also dangerous. Altimeters and ground penetratingradars (GPRs) have been used to measure the snow and ice Manuscript received June 14, 2018; revised February 17, 2019;accepted March 13, 2019. This work was supported by NASA underContract NNX15AB36G and Contract NNX17AD66G.  (Correspondingauthor: Seyedmohammad Mousavi.) S. Mousavi and K. Sarabandi are with the Radiation Laboratory, Depart-ment of Electrical Engineering and Computer Science, University of Michigan, Ann Arbor, MI 48109 USA (e-mail: mousavis@umich.edu;saraband@eecs.umich.edu).R. De Roo is with the Department of Climate and Space Sciences andEngineering, University of Michigan, Ann Arbor, MI 48109 USA (e-mail:deroo@umich.edu).A. W. England is with the College of Engineering and ComputerScience, University of Michigan, Dearborn, MI 48128, USA (e-mail:england@umich.edu).Color versions of one or more of the figures in this letter are availableonline at http://ieeexplore.ieee.org.Digital Object Identifier 10.1109/LGRS.2019.2905840 thickness [4]. However, active techniques are not well suitedto stand-off operation. In the search for a low-power remotesensing approach to measure icepack thickness, we developedwideband autocorrelation radiometry (WiBAR) [5], [6].This letter is an extension of [5], where the bare icepack thickness measurements were discussed. Here, the focus is onthe lake icepack thickness measurement using WiBAR in thepresence of a dry snow layer on top; and as such, for thepurpose of illustrating the fundamental physical principles,we will leave the issues that confound the retrieval, suchas radio frequency interference and target imperfections of absorption, and volume and surface scattering for later consid-eration in a follow-on publication. This letter is organized asfollows. In Section II, we introduce a simple forward modelingfor a lake icepack with a dry snow layer on top. We alsodiscuss the system design parameters to detect the snow andice. Section III includes the WiBAR measurements of a lakeicepack with a top dry snow layer. Section IV presents ourconclusions.II. W I BAR  FOR A  T WO -L AYER  L OW -L OSS  M EDIA WiBAR of a single layer of ice over water has beenpreviously discussed and investigated in [5]. The coherentinterference of rays traversing the slab different numbers of times gives rise to an emissivity spectrum that oscillatesaround a mean value, with local maxima at wavelengthsof constructive interference, and minima at wavelengths of destructive interference. A similar approach has been investi-gated by other researchers [7]–[10]; however, they were rarelyable to observe the quarter-wave resonances from the experi-mental results. In our WiBAR method, the emissivity spectrumis inverse Fourier transformed to obtain an autocorrelationfunction (ACF), in which the mean propagation time manifestsitself as a local maximum at a time lag greater than zero. Thedelayed ray arrives at the radiometer with the time delay  τ  p ,relative to the direct ray, as given in [5] τ  p  = 2 c    d  0   n  p (  z ) 2 −  sin 2 θ   dz  (1)where  θ   is the incidence angle,  n  p  is the refractive index of the pack,  d   is the pack’s thickness in the  z -(vertical) direction,and  c  is the speed of light in free space. The microwavepropagation time  τ  p  through the pack yields a measure of its vertical extent. 1545-598X © 2019 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.  This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 2 IEEE GEOSCIENCE AND REMOTE SENSING LETTERS Fig. 1. Remote sensing of microwave travel time through the icepack withthickness  d  ice  in the presence of a snow cover with thickness  d  snow .  A. Forward Modeling of an Icepack With a Top Snow Layer  The presence of a snowpack on the icepack adds anothermultipath. The configuration of this layered media is shownin Fig. 1. The first delayed ray of icepack and snowpack arrives at the radiometer with the time delays  τ ice  and  τ snow ,respectively, relative to the direct ray. These time delays aregiven by (1).The transfer matrix method can be used to calculate theemissivity of this medium [11], [12], as given in (2), as shownat the bottom of this page, where  ω  =  2 π  f   ,  ¯ e  is the meanemissivity over frequency,  A i ,  A s ,  A 󰀶 , and  A 󰀱  are one-half of the amplitudes of the ripple due to icepack time delay,snowpack time delay, sum of the time delays, and differenceof the time delays, respectively, in the emissivity as a functionof frequency. They are given by ¯ e  =  ( 1  − |  R 01 | 2 )( 1  − |  R 12 | 2 )( 1  − |  R 23 | 2 ) / C  0  (3a)  A i  =  R 12  R 23 ( 1  + |  R 01 | 2 ) / C  0  (3b)  A s  =  R 01  R 12 ( 1  + |  R 23 | 2 ) / C  0  (3c)  A 󰀶  =  R 01  R 23  / C  0  (3d)  A 󰀱  =  R 01  R 23 |  R 12 | 2 / C  0  (3e)where  C  0  =  1  + |  R 01 | 2 |  R 12 | 2 + |  R 01 | 2 |  R 23 | 2 + |  R 12 | 2 |  R 23 | 2 ,  R 01  is the Fresnel reflection coefficient of air to snow,  R 12  isthe Fresnel reflection coefficient of snow to ice, and  R 23  isthe Fresnel reflection coefficient of ice to water.  B. Resolving the Time Delays of the Snow and Ice Layers For two equal amplitude delay peaks in time,  τ 1  and  τ 2 ,the classic criterion for the resolution is the width of thewindow at the half power points, as explained in [5] and [13].However, since the amplitudes of the two delay peaks of iceand snow are not equal, this criteria to resolve two equalamplitude peaks is not accurate. An alternative criteria is toresolve two unequal amplitude delay peaks with an amplitudedifference,  | 󰀱  A | . If the weak amplitude delay peak residesafter the first sidelobe level of the strong amplitude delay peak,in order to resolve the two peaks, the time difference betweenthe two delay peaks  | 󰀱τ |  should satisfy | 󰀱τ |  >  max  τ main lobe 2 ,  τ FSLL  ×  2  | 󰀱  A |−| FSLL || SLF |    (4) Fig. 2. Simulated autocorrelation response with a delay peak at 1 ns.The bandwidth is 3 GHz. The delay peak has an amplitude of   − 3 dB forthe rectangular window (black solid line) and  − 15 dB for the Hammingwindow (red dashed line). where SLF (dB/Octave) is the sidelobe fall-off of the windowfunction, FSLL (dB) is the first sidelobe level of the windowfunction, and  τ FSLL  =  ( 1 / 2 )(τ main lobe + τ sidelobe )  is the locationof the FSLL peak, where  τ main lobe  =  (( 2 ζ  main lobe )/ F  s )  is themain lobe width between the first zero crossings,  ζ  main lobe  is afactor depending on the window function,  F  s  is the frequencybandwidth, and  τ sidelobe  is the width of the first sidelobe.For instance, SLF  = − 3 dB/Octave, FSLL  = − 6.5 dB, and ζ  main lobe  =  1 for the rectangular window, whereas SLF  =− 3 dB/Octave, FSLL  = − 21.5 dB, and  ζ  main lobe  =  2 for theHamming window. The width of the first sidelobe in both theHamming and the rectangular window is the same and equalto  ( 1 / F  s ) . As a case in point, to detect a delay peak 1 ns awayfrom the zero delay peak with  F  s  =  3 GHz,  | 󰀱  A |  <  9 . 5 dB,and  | 󰀱  A |  <  22 . 2 dB are required for the rectangular and theHamming window functions, respectively. This effect is shownin Fig. 2 for which the stronger delay peak is that at  τ  =  0 ns.On the other hand, if the weak amplitude delay peak residesbetween the  ( 1 / 2 )τ main lobe  and  τ FSLL  of the strong amplitudedelay peak, the two peaks can be resolved if   | 󰀱  A |  <  | FSLL | .As an example, to detect a delay peak at 0.4 ns using therectangular window,  | 󰀱  A |  should be less than 6.5 dB, whereasto detect a delay peak at 0.8 ns using the Hamming window, | 󰀱  A |  should be less than 21.5 dB, as shown in Fig. 3. In thecase of multiple delay peaks, (4) still holds for any two closelyspaced peaks in time, but the sidelobe parameters are set bythe peak with the highest sidelobes at the two close peaks,which is often the zero delay peak in the ACF.If the autocorrelation delay peaks can be resolved, theiramplitudes and positions could be biased due to the sidelobeleakage of both the zero delay and other peaks. The biasfrom the zero delay peak is visible in both Figs. 2 and 3.To reduce the effects of this bias on amplitude and time delaydetection, the window function should exhibit low-amplitudesidelobes far from the main lobe, and the transition to the low e (  f   )  =¯ e 1  +  2 [  A i  cos (ωτ ice )  +  A s  cos (ωτ snow )  +  A 󰀶  cos (ω(τ snow  +  τ ice ))  +  A 󰀱  cos (ω(τ snow  −  τ ice )) ] (2)  This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. MOUSAVI  et al. : WiBAR FOR LAKE ICEPACK THICKNESS MEASUREMENT 3 Fig. 3. Simulated autocorrelation response with delay peaks at 0.4 ns with − 2-dB amplitude for the rectangular window (black solid line) and at 0.8 nswith  − 15-dB amplitude for the Hamming window (red dashed line). Thebandwidth is 3 GHz.Fig. 4. Simulated autocorrelation response of a 35.5-cm icepack ( n ice  =√  3 . 15) without (red dashed line) and with (solid black line) a 3-cm snowpack ( ρ s  =  210 kg/m 3 ,  n snow  =  1 . 18) on top. The Hamming window was used,and the bandwidth was 3 GHz ( θ   =  0 ◦ ,  τ ice  =  4 . 2 ns,  τ snow  =  0 . 2 ns, and | 󰀱  A | ≈  1 . 4 dB). sidelobes should be very rapid [13]. On the other hand, if theautocorrelation delay peaks cannot be resolved, we can onlydetect one of the two adjacent peaks with biased amplitude andposition. The position of the detected peak would be biased tothe crossover point (halfway between the two peaks) for twoequal amplitude peaks or to a point between the crossoverpoint and the strong delay peak for two unequal amplitudepeaks, and it would get closer to the stronger peak as the  | 󰀱  A | becomes larger. For example, in the case of using a Hammingwindow with  F  s  =  3 GHz, the time delay of a 35.5-cm icepack is shifted from 4.2 to 4.3 ns in the presence of a thin 3-cmdry snowpack, as shown in Fig. 4. C. Limits of Detection The emissivity can be approximated with a Taylor expan-sion [5], for which the zeroth- and first-order terms are e (  f   )  ≈  e ( 1  −  2 [  A i  cos (ωτ ice )  +  A s  cos (ωτ snow )  +  A 󰀶 × cos (ω(τ snow + τ ice ))  +  A 󰀱  cos (ω(τ snow  −  τ ice )) ] ). (5) Fig. 5. Half amplitude of the ripples with respect to the incidence angle( ρ s  =  210 kg/m 3 ). The Brewster angle between air and snow is at about θ   =  50 ◦ , and the Brewster angle between snow and ice is at about  θ   =  80 ◦ . The absolute value of the ripple amplitudes is plotted withrespect to the incident angle shown in Fig. 5. In H-polconfiguration, it can be observed that  A 󰀱  has the lowest valueat all the angles, whereas  A i  has the largest value from nadirup to 75 ◦ . After  θ   =  75 ◦ ,  A 󰀶  becomes dominant. All theripple amplitudes in each detection scenario are lower in V-polconfiguration compared to H-pol configuration. Moreover,the amplitude values become close to zero at some anglesin V-pol configuration due to the Brewster angles between theair and snow and between the snow and ice.As in [5], the minimum  N  ind  is dependent on the receiveroperating characteristics, noise figure, and the discriminationfunction,  D (θ)  N  ind  >  2 (  Z  FA  +  Z  PD ) 2  N   f  ·  D 2 (θ)  ·  F  2 (6)where  F   is the receiver noise figure,  N   f   is the number of frequency points, and  Z  FA  and  Z  PD  are the number of standarddeviations from the magnitude of the expectation of the ACFin the absence and presence of a slab. They determine the falsealarm rate and probability of detection, respectively [5].The discrimination function depends only on the detec-tion scenario, the window function, and the incidence angle,as given by [5]  D 󰁠 (θ)  =  √  W  2 ( 0 ) W  1 ( 0 )  ( |  A 󰁠 |  e  −  SNR 󰀸 e absence )  (7)where  󰁠  =  i , s ,󰀶,󰀱 . As explained in [5],  D (θ)  must bepositive, so that the delayed peak rises above the autocor-relation noise floor. The square of the positive values of this discrimination function for lake icepack with top drysnowpack at  τ ice ,  τ snow ,  τ ice  + τ snow , and  τ ice  − τ snow  is shownin Fig. 6(a)–(d), respectively. For the snowpack, we used ρ s  =  210 kg/m 3 . It can be observed that  N  ind  needed to detectthe time delay for any of the four peaks is lower in H-polconfiguration compared to the V-pol configuration regardlessof the window function at most of the incidence angles.Moreover, similar to Fig. 5, we are only able to detect the timedelay at few angles in V-pol configuration due to the Brewsterangle between the air and snow and between the snow and  This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. 4 IEEE GEOSCIENCE AND REMOTE SENSING LETTERS Fig. 6. Square of the positive values of the discrimination function  D (θ)  for the delay peak at (a)  τ ice , (b)  τ snow , (c)  τ ice  + τ snow , and (d)  τ ice  − τ snow  as afunction of incidence angle. The snowpack density is  ρ s  =  210 kg/m 3 ( n snow  =  1 . 18). ice while we can detect the time delay at nearly all incidenceangles in H-pol regardless of the window function. The choiceof the Kaiser–Bessel ( α k   =  3 . 02, FSLL  = − 35 dB, and ζ  main lobe  =  2 . 39) window function in postprocessing couldbe helpful due to its lower sidelobes at the expense of coarserresolution. The minimum  N  ind  is inversely proportional to theripple amplitudes shown in Fig. 5, indicating that  τ ice  is theeasiest lag to measure, and  τ ice  −  τ snow  is the most difficult.Finally, Fig. 6 predicts that the microwave propagation timesin both lake icepack and the top dry snowpack are detectable atangles away from nadir to close to grazing with appropriatechoice of polarization and window function, similar to theresults shown in Fig. 5.III. F IELD  M EASUREMENTS AND  R ESULTS A snow-covered lake icepack on Lake Douglas at the Uni-versity of Michigan Biological Station (UMBS) was measuredas part of the same campaign as reported in [5]. The H-polmeasurements were made by a 7–10-GHz WiBAR instrumentwith  N  ind  =  3  ×  10 5 [5]. At these frequencies, the ice andsnow typically have negligible volume and surface scattering.The lake ice observations by WiBAR were made from about11:00  AM  until around 1:00  PM  on March 3, 2016. Theair temperature was about  − 6  ◦ C and  − 3  ◦ C at 11:00  AM and 1:00  PM , respectively. The ground truth measurementsof ice and snow thicknesses were done with a tape measureat the conclusion of the WiBAR measurements. The groundtruth values of ice and snow thicknesses are 35 . 5 and 3 . 9 cm,respectively. The snow was dry and its density was about ρ s  =  210 kg/m 3 .Among the measurements at different incidence angles, onlythe one at the largest incidence angle, 69.4 ◦ , showed distinc-tion between the two layers in the ACF. The autocorrelationresponse obtained from the measured emissivity spectra aftercalibration [5] using both rectangular and Hamming windowfunctions is shown in Fig. 7. The peaks in the ACF near4 ns are considered to be detected because they rise at leasttwo standard deviations above the expectation of the ACFin the absence of a slab. Fig. 5 shows that the ice peak isonly slightly stronger than the sum peak at this incidenceangle. The two delay peaks related to the ice thickness at τ ice  =  3 . 56 ns and the sum of the ice and the snow thick-nesses at  τ ice  +  τ snow  =  3 . 77 ns are almost resolvable usingthe rectangular window (corresponding to  d  ice  =  35 . 4 cm, d  snow  =  4 . 3 cm) while only one peak related to the icethickness is resolved at  τ ice  =  3 . 66 ns using Hammingwindow (corresponding to  d  ice  =  36 . 4 cm). There is about0.1 ns (1 cm) difference between the  τ ice  ( d  ice ) detected usingrectangular and Hamming window functions. This differenceis due to bias of the position of the detected peak, as explainedin Section II-B.  This article has been accepted for inclusion in a future issue of this journal. Content is final as presented, with the exception of pagination. MOUSAVI  et al. : WiBAR FOR LAKE ICEPACK THICKNESS MEASUREMENT 5 Fig. 7. ACF of the lake icepack with top dry snowpack measured on DouglasLake on March 03, 2016 ( θ   =  69 . 4 ◦ ,  d  ice  =  35 . 5 cm,  d  snow  =  3 . 9 cm, and ρ s  =  210 kg/m 3 ). The next terms of the Taylor expansion in (5) can alsobe observed. The delay peaks at 2 τ ice  =  7 . 24 ns and2 (τ ice  +  τ snow )  =  7 . 74 ns are detected using the Hammingwindow function, and from these, the snow thickness isinferred to be about 5.1 cm. This value is about 0.8 cm largerthan the snow thickness measured using rectangular windowfunction. This difference is due to the leakage of the sidelobes of the two close delay peaks since the ratio of the peak amplitude to sidelobe level is lower. In contrast, the rectangularwindow has such high sidelobes that these peaks cannot bedetected.IV. C ONCLUSION This letter presents the ability of the WiBAR in measuringthe snow and ice packs’ thicknesses in a lake icepack with adry snowpack on top. The presence of a snowpack can affectthe WiBAR’s capability in measuring the ice thickness bybiasing the amplitude and location of the desired delay peak if the snowpack is too thin to be resolved. This issue can bemitigated for a sufficiently thick snowpack by a proper choiceof the window function.The minimum instrument requirements to retrieve both τ ice  and  τ snow  are derived, given an ideal slab, i.e., smoothinterfaces with no absorption or volume or surface scattering.The H-polarized observations are superior to V-pol at most of the incidence angles. Detecting delay peaks is easier at  τ ice  and τ ice + τ snow  because fewer number of independent samples areneeded. Finally, detecting the ice and snow packs’ thicknessesare easier at incidence angles closer to grazing with theappropriate choice of polarization and window function.An X-band WiBAR measurement of lake ice 36 cm thick topped with a 4-cm dry snow layer retrieved propagation timesof the ice and snowpack thickness that corresponds to within2 cm of the ground truth.A CKNOWLEDGMENT The authors would like to thank R. Vande Kopple and theUniversity of Michigan Biological Station for their assistancein the field campaign.R EFERENCES[1] X. Fang and H. G. Stefan, “Long-term lake water temperature and icecover simulations/measurements,”  Cold Regions Sci. Technol. , vol. 24,no. 3, pp. 289–304, 1996.[2] P. H. Gleick,  Water in Crisis: A Guide to the World’s Fresh Water  Resources . London, U.K.: Oxford Univ. Press, 1993.[3] J. F. Manwell, C. N. Elkinton, A. L. Rogers, and J. G. McGowan,“Review of design conditions applicable to offshore wind energy systemsin the United States,”  Renew. Sustain. Energy Rev. , vol. 11, no. 2,pp. 210–234, 2007.[4] R. J. Galley, M. Trachtenberg, A. Langlois, D. G. Barber, and L. Shafai,“Observations of geophysical and dielectric properties and groundpenetrating radar signatures for discrimination of snow, sea ice andfreshwater ice thickness,”  Cold Regions Sci. Technol. , vol. 57, no. 1,pp. 29–38, 2009.[5] S. Mousavi, R. D. De Roo, K. Sarabandi, A. W. England,S. Y. E. Wong, and H. Nejati, “Lake icepack and dry snowpack thickness measurement using wideband autocorrelation radiometry,”  IEEE Trans. Geosci. Remote Sens. , vol. 56, no. 3, pp. 1637–1651,Mar. 2018.[6] A. W. England, “Wideband autocorrelation radiometric sensing of microwave travel time in snowpacks and planetary ice layers,”  IEEE Trans. Geosci. Remote Sens. , vol. 51, no. 4, pp. 2316–2326,Apr. 2013.[7] C. T. Swift, W. L. Jones, R. F. Harrington, J. C. Fedors,R. H. Couch, and B. L. Jackson, “Microwave radar and radiometricremote sensing measurements of lake ice,”  Geophys. Res. Lett. , vol. 7,no. 4, pp. 243–246, 1980.[8] C. T. Swift, D. C. Dehority, A. B. Tanner, and R. E. Mcintosh, “Passivemicrowave spectral emission from saline ice at C-band during thegrowth phase,”  IEEE Trans. Geosci. Remote Sens. , vol. GRS-24, no. 6,pp. 840–848, Nov. 1986.[9] R. F. Harrington, “The development of a stepped frequency microwaveradiometer and its application to remote sensing of the Earth,” NASA,Washington, DC, USA, Tech. Rep. TM-81847, Jun. 1980.[10] J. J. Apinis and W. H. Peake, “Passive microwave mapping of icethickness,” Electrosci. Lab., Ohio State Univ., Columbus, OH, USA,Tech. Rep. FR-3892-2, Aug. 1976.[11] L. Tsang, J. A. Kong, and R. T. Shin,  Theory of Microwave RemoteSensing . New York, NY, USA: Wiley, 1985.[12] S. Mousavi, R. De Roo, K. Sarabandi, and A. W. England, “Effectof a thin dry snow layer on the lake ice thickness measurement usingwideband autocorrelation radiometry,” in  Proc. IEEE Int. Geosci. RemoteSens. Symp. (IGARSS) , Jul. 2018, pp. 7109–7112.[13] F. J. Harris, “On the use of windows for harmonic analysis with thediscrete Fourier transform,”  Proc. IEEE  , vol. 66, no. 1, pp. 51–83,Jan. 1978.
Search
Similar documents
View more...
Tags
Related Search
We Need Your Support
Thank you for visiting our website and your interest in our free products and services. We are nonprofit website to share and download documents. To the running of this website, we need your help to support us.

Thanks to everyone for your continued support.

No, Thanks
SAVE OUR EARTH

We need your sign to support Project to invent "SMART AND CONTROLLABLE REFLECTIVE BALLOONS" to cover the Sun and Save Our Earth.

More details...

Sign Now!

We are very appreciated for your Prompt Action!

x