Home & Garden

Pan F., Peters-Lidard C.D., Sale M.J., and King A.W., A comparison of geographical information systems-based algorithms for computing the TOPMODEL topographic index. Water Resources Research, 40(6), Art. No. W06303, 2004.

[1] The performance of six geographical information systems (GIS)-based topographic index algorithms is evaluated by computing root-mean-square errors of the computed and the theoretical topographic indices of three idealized hillslopes: planar,
of 11
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
  A comparison of geographical information systems–  – basedalgorithms for computing the TOPMODEL topographic index Feifei Pan Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USA Christa D. Peters-Lidard Hydrological Sciences Branch, NASA Goddard Space Flight Center, Greenbelt, Maryland, USA Michael J. Sale and Anthony W. King Environmental Sciences Division, Oak Ridge National Laboratory, Oak Ridge, Tennessee, USAReceived 2 February 2004; revised 1 March 2004; accepted 20 April 2004; published 24 June 2004. [ 1 ] The performance of six geographical information systems (GIS)-based topographicindex algorithms is evaluated by computing root-mean-square errors of the computedand the theoretical topographic indices of three idealized hillslopes: planar, convergent,and divergent. In addition to these three idealized cases, two divergent hillslopes withvarying slopes, i.e., concave (slopes decrease from top to bottom) and convex (slopesincrease from top to bottom) are also tested. The six GIS-based topographic indexalgorithms are combinations of flow direction and slope algorithms: i.e., single flowdirection (SFD), biflow direction (BFD), and multiple flow direction (MFD) plus methodsthat determine slope values in flat areas, e.g., W-M method [ Wolock and McCabe , 1995]and tracking flow direction (TFD) method. Two combinations of horizontal resolutionand vertical resolution of the idealized terrain data are used to evaluate those methods.Among those algorithms the MFD algorithm is the most accurate followed by theBFD algorithm and the SFD algorithm. As the vertical resolution increases, the errorsin the computed topographic index for all algorithms decrease. We found that theorientation of the contour lines of planar hillslopes significantly influences the SFD’scomputed topographic index. If the contour lines are not parallel to one of eight  possible flow directions, the errors in the SFD’s computed topographic index aresignificant. If mean slope is small, TFD becomes more accurate because slope valuesin flat areas are better estimated. I   NDEX  T   ERMS  : 1899 Hydrology: General or miscellaneous;1824 Hydrology: Geomorphology (1625); 1832 Hydrology: Groundwater transport; K   EYWORDS  : GIS,TOPMODEL, topographic index, single flow direction algorithm, biflow direction algorithm, multipleflow direction algorithm Citation: Pan, F., C. D. Peters-Lidard, M. J. Sale, and A. W. King (2004), A comparison of geographical information systems–basedalgorithms for computing the TOPMODEL topographic index, Water Resour. Res. , 40 , W06303, doi:10.1029/2004WR003069. 1. Introduction [ 2 ] TOPMODEL is a topography-based concept for watershed hydrology modeling. Since the TOPMODELwas first proposed in 1979 [  Beven and Kirkby , 1979], it has been widely used to study the effects of topography onhydraulic processes including flood frequency, streamflowgeneration, flow paths, geomorphic characteristics, andwater quality [ Wolock and McCabe , 1995]. In addition tothe success of the TOPMODEL concept in traditionalhydrologic modeling, this concept has also been successfullyincorporated into several ecosystem-atmosphere modelsincluding theRegionalHydro-EcologicalSimulation System(RHESSys) [e.g., Ford et al. , 1994]; the Land Ecosystem-Atmosphere Feedback model (LEAF-2) [ Walko et al. , 2000];the TOPMODEL-based Land Atmosphere Transfer Scheme(TOPLATS) [  Famiglietti and Wood  , 1994; Peters-Lidard et al. , 1997]; the Catchment Model [  Koster et al. , 2000]; andthe Common Land Model (CLM) [  Dai et al. , 2003].[ 3 ] To apply the TOPMODEL, a modeled catchment is partitioned with a regular grid or lattice. The so-called‘topographic index’ is then calculated for each cell in thecatchment. The topographic index, ln( a /tan b ), is the naturallogarithm of the ratio of the specific flow accumulation area a to the ground surface slope tan b . The surface slope can beevaluated from digital elevation model (DEM) data. Thespecific flow accumulation area is the total flow accumula-tion area (or upslope area) A through a unit contour length  L . To compute the total flow accumulation area A , flowdirections are tracked upslope, starting from the cell of interest to the upstream divide of the watershed, and thentracked downslope accumulating cells contributing to thedrainage area of the cell of interest. Here we note anuncertainty associated with the definition of the flowaccumulation area A . From the presentations of  Beven and  Copyright 2004 by the American Geophysical Union.0043-1397/04/2004WR003069$09.00 W06303 WATER RESOURCES RESEARCH, VOL. 40, W06303, doi:10.1029/2004WR003069, 20041 of 11  Wood  [1983, Figure 2] and Kirkby [1997, Figure 1], one canreasonably conclude that the flow accumulation area isdefined along the ground surface. However, using DEMsin geographical information systems (GIS), the computedflow accumulation area is generally the area projected tox-y plane, and this calculation of  A has become standard practice. The difference between these two areas is negligi- ble if the slope is less than 0.5 (m/m), and most of the slopesin the watersheds to which the TOPMODEL is applied areless than 0.5 (m/m) [e.g., Montgomery and Dietrich , 1992].Therefore, for consistency with standard practice, we adopt the convention of calculating the flow accumulation area or upslope area as the area projected to x-y plane.[ 4 ] The slope term (tan b ) in the topographic index arisesfrom the assumption that the surface of the water table is parallel to the ground surface. Thus the local hydraulicgradient is assumed to be equal to the slope of the ground.Because flow direction depends on hydraulic gradient, or ground surface slope, flow direction and the calculation of the upslope accumulation area should be consistent with thelocal slope value that is used to compute ln( a /tan b ) . Thusthe computed topographic index is dependent upon thecalculation of both slope and flow direction.[ 5 ] Although many researchers have investigated algo-rithms for calculation of slope and flow direction, thosestudies tend to focus on either slope [e.g., Jones , 1998;  Zhangetal. ,1999]orflowaccumulationarea[e.g., Tarboton ,1997; Rieger  , 1998]. Few studies have examined the com- bined effects of slope and flow direction algorithms on thetopographic index. The works of  Quinn et al. [1991], Wolock and McCabe [1995], and Mendicino and Sole [1997] are exceptions. However, these studies only inves-tigated the difference in the statistical moments or distribu-tions of the computed topographic index; no comparisons of errors between the ‘‘true,’’ i.e., analytically solved, and thenumerically computed topographic indices were carried out.Comparing against ‘‘truth’’ is difficult because an analyticalexpression for the real terrain does not normally exist.Therefore one cannot usually determine which topographicindex algorithm is more accurate. This paper seeks toresolve that issue and provide objective evaluation of theappropriate numerical algorithm for calculating the topo-graphic index. 2. Topographic Index Algorithms [ 6 ] A topographic index algorithm is actually a combi-nation of two algorithms, one to calculate flow directionand another to calculate slope. Here we describe threecommonly used flow direction algorithms that allow flowin one, two, and more than two directions. We identify eachtopographic index algorithm by the name of the associatedflow direction algorithm (e.g., single flow direction, biflowdirection, and multiple flow direction), because each flowdirection algorithm uses its own algorithm to determineslope. 2.1. Single Flow Direction (SFD) Algorithm [ 7 ] The single flow direction (SFD) algorithm[ O’Callaghan and Mark  , 1984] calculates flow direction asthe steepest slope direction, which is determined by theMaximumDownwardGradient(MDG).ThisSFDalgorithm,also known as the D8 algorithm, is widely used in DEMdata analysis [e.g., O’Callaghan and Mark  , 1984; Band  ,1986; Greenlee , 1987; Mark  , 1988; Jenson and Domingue ,1988] and GIS software (e.g., the ‘‘FLOWDIRECTION’’function in ARC/INFO GRID).[ 8 ] MDG computes the downhill elevation gradients of a3  3 cell window along eight directions (i.e., elevation of the center cell minus elevation of each of its eight neighborsdivided by the distance between those two cells). The slopeof the central cell is calculated as the largest of the eight directions. Steepest slope direction is the direction from thecentral cell to the neighbor generating the largest downhillelevation gradient.[ 9 ] If the central cell has a lower elevation than one of itsneighbors, the downhill gradient along this direction isnegative. If the calculated slope (i.e., the largest one amongeight directions) is less than zero, this cell is called a sink or a pit. At a sink, there is only inflow, and no outflow. However,it is often computationally required to force watersheds tohave outlets. Therefore sinks in DEM data are usually filled before watersheds can be delineated and other hydrologic parameters estimated. One common filling approach for sinks generated by MDG is to raise the elevation of a sink to the lowest elevation among its neighbors [e.g., Jenson and  Domingue , 1988; Tarboton et al. , 1991]. This algorithm isused in the ARC/INFO GRID (i.e., the ‘‘FILL’’ function).For example, the recalculated slope value of the ‘‘filled’’ cellis now zero; the cell and its neighbors form a flat area. Todetermine flow direction in a flat area, the method suggested by Jenson and Domingue [1988] (e.g., the ARC/INFOFLOWDIRECTION function) is used.[ 10 ] Once the single flow direction is determined for eachcell, flow accumulation area (or upslope area) (  A ) iscalculated using a recursive procedure [e.g., Tarboton ,1997]. The specific flow accumulation area ( a ) is (  A )divided by a contour length, which is equal to the grid sizeor horizontal resolution of the DEM. The slope is set to bethe maximum downward elevation gradient. 2.2. Biflow Direction (BFD) Algorithm [ 11 ] The biflow direction algorithm (BFD), also known asD 1 , was proposed by Tarboton [1997]. In this algorithmthe 3  3 cell window is divided into 8 triangular facets.The slope direction and magnitude of each facet arecompared. The steepest downward direction is chosen anddivided into two directions along the edges forming that facet. The proportion of flow along each edge is inversely proportional to the angle between the steepest downwarddirection and the edge. Therefore at most two flow direc-tions can be assigned to each cell. The contour length isdefined as the grid cell size (DEM resolution), and the slopeis set to be the largest slope of 8 facets [ Tarboton , 1997].[ 12 ] Like the SFD algorithm, the BFD algorithmhas trouble defining flow direction in flat areas.A GIS-based algorithm for returning flow directions[  Jenson and Domingue , 1988], such as the ARC/INFO‘‘FLOWDIRECTION’’ function, can be used to determinethe flow direction in flat areas, such that only one flowdirection is assigned to each cell in a flat area. 2.3. Multiple Flow Direction (MFD) Algorithm [ 13 ] Quinn et al. [1991] first suggested the multiple flowdirection (MFD) algorithm to improve representation of theconvergence or divergence of flow. Wolock and McCabe 2 of 11 W06303 PAN ET AL.: A COMPARISON OF TOPOGRAPHIC INDEX ALGORITHMS W06303  [1995] showed how to implement this algorithm usingARC/INFO GRID functions. Unlike the SFD and BFDalgorithms, the MFD algorithm allows flow in more thanone or two directions. In the MFD algorithm, each flowdirection is weighted by the downward elevation gradient (i.e., from the central cell to each of its 8 neighbors)multiplied by a ‘‘contour length’’. There are two ways toset the contour length: i.e., (HR/2) and (  ffiffiffi 2 p  HR /4) [ Quinnet al. , 1991], or 0.6HR and 0.4HR [ Wolock and  McCabe , 1995], for cardinal and diagonal flow directions,respectively, where HR is the horizontal resolution, 0 ° gridcell size of the DEM. Here, we assume the HR in eastingand northing directions are the same (which is true for most DEM data). To differentiate those two algorithms for assigning the contour lengths, we call Quinn’s method theMFD algorithm, and Wolock and McCabe’s algorithm theMFD* algorithm. The slope is the weighted summation of all positive downward gradients [ Quinn et al. , 1991]. As for the MFD algorithm, a GIS-based algorithm for returningflow directions [  Jenson and Domingue , 1988], such asARC/INFO flow direction function can used to determinethe flow directions in flat areas. 2.4. Recalculation of Slopes in Flat Areas [ 14 ] Because the topographic index cannot be calculateddirectly in flat areas, estimating slopes in flat areas is acritical issue. Wolock and McCabe [1995] defined theminimum slope value to be equal to (0.5  verticalresolution)/(horizontal resolution). We refer to this as theW-M method.[ 15 ] To maintain consistency between flow direction andslope in flat areas, we propose a simple algorithm called thetracking flow direction (TFD) method. This algorithm startsfrom cell X  and following the flow direction determined byARC/INFO ‘‘FLOWDIRECTION’’ function that searchesfor the nearest cell Y  with a lower elevation than cell X  . Theslope at cell X  is then set to be the ratio of elevationdifference to length of flow path between cells X  and Y  (see Figure 1). If a nonzero slope value cannot be found for the cell over the whole domain, the value of  À 1 is assignedto this cell, indicating that slope is undefined. The minimumof all defined slopes is then assigned to each undefined cell. 3. Idealized Hillslopes [ 16 ] To compare the computed and theoretical ‘‘true’’topographic indices, we apply six GIS-based topographicindex algorithms to three idealized hillslopes, i.e., planar,divergent and convergent, for which we can calculate thetheoretical topographic index analytically. Although thosehillslopes are ideal simplifications, they represent threetypes of real terrain: terrace, divergent and convergent.[ 17 ] For the planar and convergent hillslopes, we use oneslope value for each case, i.e., slope is spatially uniform. For the divergent hillslopes, in addition to the constant slopecases, we also test two divergent hillslopes with varyingslopes, i.e., concave (slope decreases from top to bottom),and convex (slope increases from top to bottom). Althoughslopes are not spatially uniform for those concave andconvex divergent hillslopes, to simplify the problem, weset slope constant along each contour line. Under suchconditions, we can theoretically compute the specific flowaccumulation area of each cell along a contour line as thetotal upslope area divided by the length of the contour line.This calculation is valid only when the contour line isclosed for the divergent or convergent hillslopes, or thecontour line is infinitely long for the planar hillslopes.[ 18 ] Prior to introducing each idealized case, we first define three domains. The computation domain is definedas the region where we use a recursive algorithm to computethe flow accumulation area based on flow directions deter-mined by the SFD, or BFD or MFD algorithms. We canthink of the computational domain as a part of a watershed.However, this part is not just an arbitrary domain. It startsfrom a divide or a local peak and covers a part of thedownstream area, since we will set the top boundary parallelto the contour lines. Therefore, when we compute the flowaccumulation area, we do not miss any contribution fromupslope crossing the top boundary.[ 19 ] The influence domain is defined as the area insidethe computational domain that receives the flow contribu-tion from the lateral boundaries. The comparison domain isinside the computational domain but with the influencedomain excluded. Our theoretical calculation of the topo-graphic index is valid only if there is no influence from thelateral boundaries. Therefore, to compare the computed andthe theoretical topographic indices properly, we need tocarry out such comparison inside the comparison domain, because only inside the comparison domain are the com- puted flow accumulation area and thus the computedtopographic index not influenced by the lateral boundaries. 3.1. Planar Hillslopes [ 20 ] A three-dimensional, idealized planar hillslope isillustrated in Figure 2. In a NROW   NCOL domain(  NROW  and NCOL are numbers of rows and columns),the coordinates of each cell are given by:  x ¼ j  À 0 : 5 NCOL À 1 ð Þ½  * HR ;  j  ¼ 0 ; . . . ;  NCOL À 1  y ¼ 0 : 5* NROW  À 1 ð ÞÀ i ½  * HR ; i ¼ 0 ; . . . ;  NROW  À 1 ð 1 Þ where (i, j) are row and column indices of each cell, HR isgrid cell size or horizontal resolution of DEM data. In thisstudy, we set  NROW  = 2001, and NCOL = 4001. Figure 1. An example showing the TFD method for recalculating slope values in flat areas. Numbers in cells areelevation values. Following the ARC/INFO derived flowdirections, we identify the closest cell to cell X with a lower elevation is cell Y. The slope value for cell X is computed asthe ratio of the elevation difference to the length of flow path between X and Y, which is 2/(2  ffiffiffi 2 p  HR + HR ), where  HR is the cell size. W06303 PAN ET AL.: A COMPARISON OF TOPOGRAPHIC INDEX ALGORITHMS3 of 11 W06303  [ 21 ] The theoretical elevation Z  t  of each cell is given by:  Z  t  ¼ S  *  slope ¼ S  *tan b ð 2 Þ where tan b is slope, b is slope angle (0 ° < b < 90 ° ), and S  isthe distance from (x o , y o ) to the contour line passing cell(x,y) given by S  ¼ x À  x o ð Þ tan a þ y À  y o ð Þj j  ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tan a ð Þ 2 þ 1 q  ;  x o ¼À 0 : 5 NCOL À 1 ð Þ  HR ;  y o ¼À 0 : 5 NROW  À 1 ð Þ  HR ; ð 3 Þ where a is the angle between each contour line and x axis(from the contour line to x axis in the counterclockwisedirection) and we call it contour angle (Figure 2).[ 22 ] To investigate the effects of the vertical resolution(VR) of elevation data on the errors in the computedtopographic index, we set elevation z (which is used for determining flow directions and slopes) as follows:  z  ¼ VR *int  Z  t  ð Þ ; if fractional part of  Z  t  < 0 : 5 VR * 1 þ int  Z  t  ð Þ½  ; if fractional part of  Z  t  ! 0 : 5 & ð 4 Þ whereint()isafunctiontoconvertafloatingpointvaluetoaninteger by cutting off all fractional part. Using expressionsshown in equation (4), we can mimic the vertical resolutionof real DEM data. For example, a theoretical elevationis 125.253m. According to equation (2), the elevation is125 m at 1-m VR, or 125.3 m at 0.1-m VR, or 125.25 m at 0.01-m VR.[ 23 ] The computational domain for the planar hillslopes isconfined by four boundaries (Figure 2), i.e., TB, BB, LBand RB. TB is called the top boundary and BB is the bottom boundary. They are parallel to the contour lines and asso-ciated with the highest and lowest elevations in the com- putational domain, respectively. LB and RB are two lateral boundaries. Facing the direction of downward gradient of elevation (shown as a block arrow in Figure 2), LB is on theleft side, and RB is on the right side. We set the computa-tional domain for the planar hillslopes as follows: d  1  0 : 1 NCOL À 1 ð Þ  HR ; d  2  0 : 1 NCOL À 1 ð Þ  HR ; d  3 < 0 : 1 NCOL À 1 ð Þ  HR ; d  4  0 : 1 NCOL À 1 ð Þ  HR ð 5 Þ where d  1 , d  2 , d  3 , and d  4 are distances from srcin (0, 0) tofour boundaries: LB, RB, TB and BB.[ 24 ] In the computational domain, the theoretical topo-graphic index is given by: ln a = tan b ð Þ¼ lnupslope area = contour lengthslope ! ¼ ln D *  L L tan b ! ¼ ln D tan b ! ð 6 Þ where D is the distance from cell (x, y) to TB, and L is thecontour length. As we mentioned in section 1, the flowaccumulation area of the upslope area is the area projectedto x-y plane, i.e., DL . Although L does not appear in thefinal expression, equation (6) is valid only as L is infinite, or there is no influence due to the lateral boundaries. However,  L cannot be infinite. Therefore a proper comparison between the computed and theoretical topographic indexcan only be carried out in a comparison domain where thelateral boundaries have no influence.[ 25 ] Our first step is to remove the influence domain fromthe computational domain. In order to determine the influ-ence domain, we need to follow flow directions flowing intothe computational domain from two corner cells on the top boundary (TB). The MFD algorithm has a larger influenceregion than the SFD algorithm or the BFD algorithm, because it allows more flow dispersion. Therefore wechoose the MFD’s comparison domain as the commoncomparison domain, because it is in both SFD’s and BFD’scomparison domains.[ 26 ] Figure 3a shows 8 possible flow directions at a cell.The angles between these flow directions and x axis are 0 ° (direction 1), 45 ° (direction 2), . . . , and 315 ° (direction 8).Figures 3b, 3c, and 3d illustrate the computational, influ-ence and comparison domains for three cases correspondingto three different contour angles, i.e., 0 ° , 45 ° and 22.5 ° . Thearea between LB and AI is the influence domain due to thelateral boundary LB. The angle formed by LB and AI (i.e., aa ) is the maximum angle among the angles formed by LBand each flow direction at point A which points to thecomputational domain. Similarly, we can determine the Figure 2. Top and side views of a planar hillslope. Thecomputational domain is formed by four boundaries (LB,RB, TB, and BB). Inside the computational domain, theunshaded area is the comparison domain, and the shadedarea is affected by the lateral boundaries (LB and RB). The block arrow is the normal direction of the contour lines andalong the downward elevation gradient direction. The angle between each contour line and x axis is a , S is the distancefrom point (x o , y o ) to a contour line, D is the distance between a cell and the top boundary TB, and b is slopeangle. 4 of 11 W06303 PAN ET AL.: A COMPARISON OF TOPOGRAPHIC INDEX ALGORITHMS W06303  influence domain due to the lateral boundary RB, which isformed by RB and BI. The angle between RB and BI is bb .In Figure 3b, both points A and B have three flowdirections, i.e., 6, 7, and 8. However, only direction 8 at  point A and direction 6 at B point into the computationaldomain. Therefore angles aa and bb are equal to 45 ° .Similarly, in Figure 3c, both A and B have three flowdirections, i.e., 5, 6, and 7. Among those directions, onlydirection 7 at point A and direction 5 at point B flow intothe computational domain. Therefore angles aa and bb areequal to 45 ° . In Figure 3d, both A and B have four flowdirections, i.e., 5, 6, 7 and 8. Among those directions, onlydirections 7 and 8 at A, and directions 5 and 6 at B flow intothe computational domain. Comparing the angles, we canfind that angle aa is formed by LB and direction 8, andangle bb is formed by RB and direction 5. Therefore aa and bb are equal to 67.5 ° . Once aa and bb are determined, wecan determine the influence domain and thus the compar-ison domain. 3.2. Divergent Hillslopes3.2.1. Divergent Hillslopes With Constant Slopes [ 27 ] A top and a side view of a divergent hillslope areillustrated in Figure 4. In a NCOL  NROW  domain, thecoordinates of each cell are the same as the planar hillslopes (i.e., equation (1)). The theoretical elevation  Z  t  is given by:  Z  t  ¼  f   1 *  HR À r  *  slope ¼  f   1 *  HR À r  *tan b ð 7 Þ where tan b is slope, b is slope angle (0 ° < b < 90 ° ), r  =  ffiffiffiffiffiffiffiffiffiffiffiffiffiffi  x 2 þ  y 2 p  is radius of the contour line passing (  x , y ), and  f   1*  HR is the elevation at the center of the domain. In thisstudy, we set  NCOL = 4001, NROW  = 2001, f   1 = 3000, andthe computational domain inside the circle (exclude thesrcin (0, 0)) with a radius equal to R (Figure 4), whichsatisfies the following conditions: 0 < r    R ð 8 Þ where we set  R = 500HR in this paper.[ 28 ] In the computational domain, the theoretical topo-graphic index is given by: ln a = tan b ð Þ¼ lnupslope area = contour lengthslope ! ¼ ln p r  2 2 p r  tan b ! ¼ ln r  2tan b ! ð 9 Þ The comparison domain is the same as the computationaldomain, since the contour lines are closed. 3.2.2. Divergent Hillslopes With Varying Slopes [ 29 ] Figure 5 illustrates two divergent hillslopes withvarying slopes: slope decreases from to top bottom (con-cave, Figure 5a) and slope increases from top to bottom(convex, Figure 5b). The theoretical elevation Z  t  is given by Concave Z  t  ¼ c 1 r  À c 2 ð Þ 2 Convex : Z  t  ¼ c 3 À r  = c 4 ð Þ 2 ð 10 Þ Figure 3. (a) The 8 flow directions, i.e., 1, 2, . . . , 8. (b) Computational, influence, and comparisondomains as the angle between each contour line and x axis is 0 ° . The computational domain is formed bythe top boundary TB, bottom boundary BB, and two lateral boundaries LB and RB. The influencedomain is the shaded area formed by LB, BB, RB, BI, and AI. The comparison domain is the unshadedarea inside the computational domain with the influence domain excluded. (c) Same as Figure 3b, except the contour angle is 45 ° . (d) Same as Figure 3b, except the contour angle is 22.5 ° . W06303 PAN ET AL.: A COMPARISON OF TOPOGRAPHIC INDEX ALGORITHMS5 of 11 W06303
Similar documents
View more...
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

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!