Next Article in Journal
Morphogenesis of a Floodplain as a Criterion for Assessing the Susceptibility to Water Pollution in an Agriculturally Rich Valley of a Lowland River
Next Article in Special Issue
Glacial Lake Detection from GaoFen-2 Multispectral Imagery Using an Integrated Nonlocal Active Contour Approach: A Case Study of the Altai Mountains, Northern Xinjiang Province
Previous Article in Journal
Ensemble Kalman Filter Assimilation of ERT Data for Numerical Modeling of Seawater Intrusion in a Laboratory Experiment
Previous Article in Special Issue
Land Water-Storage Variability over West Africa: Inferences from Space-Borne Sensors
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Unstructured-Mesh Terrain Analysis and Incident Solar Radiation for Continuous Hydrologic Modeling in Mountain Watersheds

1
Department of Geography and Environmental Sustainability, University of Oklahoma, Norman, OK 73019, USA
2
Department of Civil and Architectural Engineering, University of Wyoming, Laramie, WY 82071, USA
3
Center for Automated Sensing and Sampling, School of Meteorology, University of Oklahoma, Norman, OK 73019, USA
*
Author to whom correspondence should be addressed.
Current address: 100 E. Boyd St, Norman, OK 73019, USA
Water 2018, 10(4), 398; https://doi.org/10.3390/w10040398
Submission received: 26 February 2018 / Revised: 20 March 2018 / Accepted: 22 March 2018 / Published: 28 March 2018
(This article belongs to the Special Issue Applications of Remote Sensing and GIS in Hydrology)

Abstract

:
This article presents a methodology for estimating total incoming solar radiation from Triangular Irregular Network (TIN) topographic meshes. The algorithm also computes terrain slope degree and aspect (slope orientation) and accounts for self shading and cast shadows, sky view fractions for diffuse radiation, remote albedo and atmospheric backscattering, by using a vectorial approach within a topocentric coordinate system establishing geometric relations between groups of TIN elements and the sun position. A normal vector to the surface of each TIN element describes its slope and aspect while spherical trigonometry allows computing a unit vector defining the position of the sun at each hour and day of the year. Sky view fraction, useful to determine diffuse and backscattered radiation, is computed for each TIN element at prescribed azimuth intervals targeting the steepest elevation gradient. A comparison between the sun zenith angle and the steepest gradient allows deciding whether or not the pivot element is shaded. Finally, remote albedo is computed from the sky view fraction complementary functions for observed albedo values of the surrounding terrain. The sensitivity of the different radiative components to seasonal changes in atmospheric transmissivitties and surrounding albedo is tested in a mountainous watershed in Wyoming. This methodology represents an improvement on the current algorithms to compute terrain and radiation values on unstructured-mesh terrain models. All terrain-related features (e.g., slope, aspect, sky view fraction) can be pre-computed and stored for easy access into a subsequent, progressive-in-time, numerical simulation.

1. Introduction

Terrain properties such as slope gradient and aspect are fundamental in controlling a number of point-surface, and shallow subsurface, hydrologic processes including runoff rates and generation mechanisms [1,2], infiltration and wetting front development [3,4], weathering and erosion rates [5,6] and incident solar radiation [7], among other ecohydrologic processes. On a clear day, the solar radiation flux striking an inclined surface consists of direct, diffuse, backscattered and landscape reflected (albedo) radiation [3,8]. In flat terrain, the first three components dominate. Nevertheless, in rough topography, remote shading and albedo from surrounding slopes become significant [9,10]. The total incident solar radiation and its variability in topographically complex areas is responsible for controlling magnitude and spatial variability of turbulent heat fluxes, snow ablation rates, soil moisture losses and skin temperature, among other physical variables. The uneven distribution of snow and vegetation on contrasting mountain flanks is a clear evidence of this topographically-induced variability [11,12,13].
TINs have become relevant Digital Elevation Models (DEMs) for their versatility to represent, store, access and update terrain data. In some cases, Voronoi polygons [14,15] are derived from TINs through the grouping of neighboring triangles sharing common triangle vertices (i.e., central nodes). Each voronoi element constitutes the smallest computational unit with its own slope degree and aspect determined by the steepest-gradient triangle edge downstream from each central node [16,17,18]. This procedure, however, lumps the polygon inter-triangle terrain variability into a single representative value. More traditional grid-based DEMs develop their computational algorithms on squared elements that use a 3 × 3 kernel or moving window displaced along the grid cells to compute slope and aspect [19,20,21,22,23,24,25]. This approach has been shown to underestimate slope [26,27] in highly irregular terrains, making it less appealing for hydrologic and geomorpologic studies that deal with runoff generation, flash-floods and landslides prediction. Similarly, slope aspect gets constrained to this zonal (nxn cells) algorithm conditioning its occurrence to (n × n) − 1 possible occurrences on the X-Y coordinate system. Corripio [28] attempted to resolve this for grid-based model domains by subdividing each DEM pixel in two triangles whose normal vector defines its direction in the X-Y cartesian coordinate system. However, having twice the number of domain elements, in a regularly-spaced array, increases computational burden.
A realistic representation and estimation of terrain properties and incident solar radiation in unstructured (i.e., non-grid) based model representations is primal for an accurate understanding and prediction of water and energy balances in both large-scale and hyper-resolution land/atmosphere interaction models [29]. Computing incident solar radiation in unstructured meshes has significantly advanced, in the last decade. Montero et al. [30] use an adaptive procedure of domain refinement/ derefinement based on nested meshes to compute shading. This approach, however, may become impractical for continuous hydrologic modeling as it re-creates domain meshes to facilitate the calculation of radiation and shading. The work by Marsh et al. [9], based on Clarke [31] and Montero et al. [30], relies on an Euler rotational algorithm applied on each TIN element to determine remote shading in a region of the Canadian Rockies. Nonetheless, the algorithm is applied to each single element without horizon visibility constraints to help reduce the search time. Other advances in the estimation of incident radiation for TIN domains in hydrologic models are still based on ancillary DEMs to implement a 8-direction algorithm across the entire grid domain [10] incurring extensive calculations. Alternatively, the implementation of coarse resolution radiation products (e.g., 4 km pixel resolution) from weather forecasting models, satellite imagery or radiometer-interpolated maps represents a gain in spatial cover at the expense of reduced resolution and variability [32,33,34]. There is a clear need for unstructured TIN-based methods to compute surface terrain properties and energy inputs, to take full advantage of this improved geometric representation.
In this article, we describe the development and testing of a TIN-based algorithm to compute terrain parameters and clear-sky incident solar radiation for implementation in continuous and distributed hydrologic models that employ unstructured TIN geometry. The algorithm was developed as a module for the ADHydro model [35], although it could also be partially or totally adapted by other TIN-based models. It accounts for direct solar with self-shading and cast shadows, sky-view diffuse, atmospheric backscattered and surrounding albedo-reflected components without refining or derefining the original mesh. The algorithm is presented as a standalone code whose input is a set of triangle vertex points and their elevation, surface albedo, air temperature and relative humidity. To develop and test our algorithm, we selected the Green River headwater catchment, a topographically-complex watershed in west-central Wyoming. Use of this estimation strategy can result in benefits for improved hydrologic predictions, model parsimony and reduced computational loads for serial and parallel computing in irregular-mesh terrain-based models.

2. TIN Domain Geometry and Application Catchment

The upper Green River basin in Wyoming, U.S, chief tributary to the Colorado River, is selected to conduct this study (Figure 1). It rises in western Wyoming in the Bridger-Teton National Forest and flows south through Sublette County and western Wyoming. The selected catchment has a drainage area of 1220 km 2 with a mean slope of 13.24%. The east portion of the catchment is mountainous and characterized by high elevations whose snowmelt supports the majority of water-year streamflow [36]. The presence of a large elevation relief, steep slopes and wide-changing valleys makes this an appealing region to apply and test this methodology.
The algorithm presented in this manuscript is based on an irregularly-spaced array of triangles on which a series of vectorial calculations are made. For this manuscript, a TIN was generated from a 876,215 cell, 30-m resolution DEM from the Shuttle Radar Topography Mission (SRTM) at 1 arcsec resolution (∼30 m). No further error analyses were conducted to assess the accuracy of this product. The TIN used uniform random sampling to preserve 5% the original number of nodes. The final number of TIN computational elements was not determined based on the speed or accuracy of the model results presented in this study. Nonetheless, sampling to reduce the number of original elements in 95% will represent a significant reduction in the computational burden of any hydrologic model [37,38]. The Figure 2A,B illustrate the used TIN representation.
As illustrated by an example case in Figure 2C, the smallest computational surface of a TIN is the triangle depicted by coordinates (x 1 , y 1 , z 1 ), (x 2 , y 2 , z 2 ) and (x 3 , y 3 , z 3 ). The triangle defined by these three points in space determines a normal vector by computing the cross product between the vectors u and v as shown in Equation (1). Selection of vector pairs u and v is such that the right-hand rule results in a k-positive, upward-looking, vector.
n = u × v = ı ^ j ^ k ^ x 2 x 1 y 2 y 1 z 2 z 1 x 3 x 2 y 3 y 2 z 3 z 2
The solution of this determinant provides the components of the normal vector to the example triangle in Figure 2C in terms of coordinate and elevation differences, as illustrated by Equation (2).
n = [ z 3 ( y 2 y 1 ) + z 2 ( y 1 y 3 ) + z 1 ( y 3 y 2 ) ] ı ^ [ z 3 ( x 1 x 2 ) + z 2 ( x 3 x 1 ) + z 1 ( x 2 x 3 ) ] j ^ [ z 3 ( x 2 x 1 ) + y 2 ( x 1 x 3 ) + y 1 ( x 3 x 2 ) ] k ^
Thus, dividing the components of n by its magnitude, a normal unit vector ( n u ) is obtained. Figure 3A,B illustrate the horizontal components of n u ( n u x and n u y ) for the Green River Basin along with quartile metrics of the spatially-distributed values. It is observed that west and south facing slopes slightly dominate across the landscape of this basin with marked direction contrasts on the eastern Wind River Range mountains (Figure 3A,B).

3. Data and Methods

3.1. Slope Degree and Aspect Calculation

Slope degree ( ζ ) and aspect ( α ) are derived from Equations (3) and (4):
ζ = c o s 1 n u z
α = n u x | n u x | π 2 t a n 1 n u y n u x
where the slope ( ζ ) is the angle formed between the steepest-gradient tangent to each triangle’s surface and the horizon plane such that 0 ζ 90 and aspect ( α ) is the compass direction of each slope relative to the north direction (i.e., Azimuth) such that 180 α 180 . Spatial distribution of ζ and α is illustrated in Figure 4A,B.

3.2. Clear-Sky Components of Incident Radiation to a Non-Vegetated Surface

The instantaneous, clear-sky incident radiation on an inclined surface (K i n ) is estimated as the sum of the direct beam (K d i r ), diffuse (K d i f ), atmospheric backscattered (K b a c ) and reflected by surrounding albedo (K a l b ) components [3,10] as illustrated in Equation (5).
K i n = K d i r + K d i f + K b a c + K a l b
Each of these paramters are related to a different component of the total downwelling shortwave radiation as explained in Section 3.3, Section 3.4, Section 3.5 and Section 3.6.

3.3. Direct Radiation, K d i r

The direct radiation component to an inclined TIN element depends on the cosine of the angle ( θ ) between the normal unit vector (n u ) to the surface and the unit solar vector(s), the capacity of the atmosphere to transmit the solar beam and the shading produced by prominent, surrounding topography.

3.3.1. Solar Unit Vector and Self Shading

A vector defining the sun position is obtained from the solar declination, Julian date and hour angle for each triangular element using spherical trigonometry [3,28,39]. Different solar vectors per TIN element and time step account for watershed’s geographic differences and their effect on the apparent location of the sun, a capability that is particularly useful in large-size basins. The details on the solar vector calculation are presented in Appendix A. Using this approach a unit solar vector s is computed by dividing each component by its vector magnitude. The third component of the sun vector, s z , is positive when sun is above the horizon and negative otherwise. The direct radiation component to an inclined TIN element is proportional to the cosine of the angle ( θ ) between the normal unit vector ( n u ) of the surface and the unit solar vector ( s ), as dictated by the Lambert cosine law [40,41] and expressed as the dot product between two vectors in Equation (6).
c o s ( θ ) = s · n u
The sign of this dot product determines if the sun is visible from each element ( | θ | ≤ 90) or if element self-shading occurs ( | θ | ≥ 90).

3.3.2. Shadow Casting

Shadow casts play a fundamental role in complex topography particularly during early or late day hours when the sun vector is shadowed by surrounding mountains. In order to check for remote shading to each triangle element, a shading search algorithm was developed according to an Earth curvature visibility threshold (i.e., d v ) as expressed in Equation (7):
d v = h 1 2 h 0 2 + 2 R ( h 1 h 0 )
where h 1 and h 0 are the elevations of the shading and shaded nodes, respectively and R is the Earth’s mean radius. The algorithm scans for highly visible, sky-blocking triangle vertices using an azimuth interval-search method from each triangle center T 0 , that serves as pivot, to horizon-visible (i.e., d < d v ), triangle vertices (including each pivot’s own triangle vertices; see Figure 5A).
The method divides the horizon search domain into n slice directions (e.g., u = 1, 2, 3, 4,…, n = 64) per pivot point. One blocking vertex per slice u is determined as indicated by the maximum positive elevation gradient (i.e., elevation difference-to-distance ratio) Δz/d, where Δz and d are the elevation difference and horizontal distance between the pivot and each possible sky-blocking element (Figure 5B). This ratio ensures the comparison of far and tall geomorphologic features to closer, shorter topographic objects. Having a single blocking element per slice u, the zenith angle (f v s u ) of a sight line between the pivot T 0 and the blocking TIN vertex (( Δ Z/d) m a x ) is calculated per search slice as in Equation (8) (See Figure 5C).
f v s u = a t a n d Δ Z
For each time step, the solar vector provides its horizontal azimuth ( α s ) and zenith angle ( ϕ δ ) as described by Equations (9) and (10)
α s = s x | s x | π 2 t a n 1 s y s x
ϕ δ = a t a n s x 2 + s y 2 s z 2
Thus, for a particular pivot, within the slice that contains α s , if f v s u < ( ϕ δ ) then this pivot is remotely shaded by the blocking TIN vertex at this particular time step.

3.3.3. Atmospheric Solar Beam Attenuation

Under clear-sky conditions, before reaching the ground, the total direct incident radiation suffers atmospheric scattering and absorption by air components, water vapor, dust and solid aerosols [3,42]. These effects are accounted for using transmissivity values that depend on the precipitable water content and optical air mass. The reader is referred to Appendix A for a detailed description of atmospheric transmissivity calculations. Transmissivity values were computed from hourly records of air temperature and relative humidity taken from the Upper Green River Valley Weather Station (KWYCORA2) with coordinates 110.02 W, 43.22 N, elev = 2331 m, located inside the watershed divide, nearby the City of Cora, Wyoming. On an sloping surface the total direct radiation (K d i r ) can be computed as Equation (11).
K d i r = k s D s ( t ) c o s [ θ ( t ) ] τ w a ( t ) τ d a ( t ) τ w s ( t ) τ r s ( t ) τ d s ( t ) , s z > 0 K d i r = 0 , s z 0
where k s is the solar constant (k s = 1361 w/m 2 ), D s (t) is the sun-earth distance function (0 ≤ D s ≤ 1), s z is the vertical component of the unit solar vector at each TIN element and τ w a ( t ) , τ d a ( t ) , τ w s ( t ) , τ r s ( t ) , τ d s ( t ) are the transmissivity functions resulting from the absorptive properties of water vapor, dust and solid aerosols and light scattering due to water vapor, air molecules (Rayleigh Scattering) and dust and solid aerosols, respectively [3] (see Appendix A).

3.4. Diffuse Radiation, K d i f

The fraction of the total incident radiation that is scattered from the direct solar beam by molecules or particles in the atmosphere is the diffuse radiation (k d i f ). The flux of diffuse radiation to each TIN element depends on the sky view fraction f v s . Having f v s , k d i f can be computed in a similar way as suggested by Dingman [3]. The diffuse radiation only exists when the sun is above the horizon and it is zero otherwise. This assumption neglects the effects of remaining diffuse radiation after the sunset or preceding the sunrise. The proportionality factor C was found to be equal to 0.05.
K d i f = C k s D s ( t ) τ w a ( t ) τ d a ( t ) ( 1 τ w s ( t ) τ r s ( t ) τ d s ( t ) ) f v s , s z > 0 K d i f = 0 , s z 0

3.4.1. Sky View Fraction

The visible portion of the sky from any point on the ground is called the sky view fraction (f v s ). An ideally flat landscape with no visible mountains surrounding would have f v s = 1 while a completely obstructed condition would lead to f v s = 0. A triangle element on a tall mountain summit could have f v s > 1, meaning that more than one hemisphere could be visible from that location. f v s is an important parameter for the estimation of the diffuse (K d i f ) and atmospheric backscattered (K b a c ) radiation components [8,28,43]. The complementary value f v l , also called land view fraction, indicates how much ground can be seen from a particular element. f v l = 1 − f v s f v s ≤ 1. f v l is particularly useful when calculating the total energy flux reflected by surrounding areas with high albedo (K a l b ), a significant input of energy in narrow valleys under snow cover conditions [44,45,46].
Typically, f v s is computed as the zenith angle of the horizon integrated across an n-direction search method (n = 8 for squared elements) and then averaged for all the compass directions to compute the ratio of visible to blocked sky hemisphere [7,47]. We propose to use the strategy developed in Section 3.3.2 (i.e., shadow casting) to integrate the fraction of the sky viewed from each pivot T 0 . Computationally, the calculations would not represent an additional load since an n-directional matrix of sky-blocking elements was previously created and stored. Thus, a single f v s is obtained by averaging all f v s u for the n slice search directions and dividing by 90 of full half-horizon visibility along each direction. A spatially-distributed map of f v s is illustrated in Figure 6. Note that narrow-valley floors in areas surrounded by tall and close hillslopes hold relatively low f v s values, while both flat and high elevation regions hold f v s values near unity.

3.5. Atmospheric Backscattered Radiation, K b a c

The fraction of the short wave radiation that is re-reflected back from the atmosphere (K b a c ) as a result of the land surface albedo ( α ), is computed in an analogous way to the diffuse component as expressed by Equation (12). Basin-average, seasonal albedo values were extracted from the NLDAS model [48] at a resolution of 12.5 km. This expression, however, neglects higher-order reflection terms due to successive reflections between the land-surface and the atmosphere.
K b a c = α ( K d i r + K d i f ) C τ w a ( t ) τ d a ( t ) ( 1 τ w s ( t ) τ r s ( t ) τ d s ( t ) ) f v s

3.6. Landscape Albedo Reflected Radiation, K a l

The fraction of the short wave radiation that is reflected back from the surrounding landscape’s albedo (K a l b ), can be computed as described by Equation (14).
K a l b = α K d i r ¯ f v l
where α K d i r ¯ represents the mean albedo-reflected energy from the surrounding landscape to each pivot element and f v l is the land view fraction. This approach constitutes a simplified, economic version for use in large watershed modeling applications. However, it ignores higher-order terms due to the successive reflection of energy between landscape elements and subsequent backscattered energy from the atmosphere.

4. Results

4.1. Instantaneous Flux of Incident Solar Radiation

The resulting maps of instantaneous solar radiation flux (K i n ; w/m 2 ), after including all the components from Equation (5) at each TIN element are illustrated in Figure 7A,B for two different hours (12 PM and 5 PM LST) on 20 March (Spring equinox) of year 2016. Contrasting solar radiation fluxes are observed at the two hours (12 PM and 5 PM) with the 12 PM condition (Figure 7A) receiving significantly more radiation across the entire basin, particularly for south-facing elements. At 5 PM (Figure 7B), west-facing slopes are the only landscape features that maintain comparable radiation values to the 12 PM condition. In both cases, the high variability of radiation fluxes over the east portion of the watershed reveals the shading effects by nearby, steep slopes and predominance of significant topographic heterogeneity.

4.2. Slope Degree and Aspect Controls on Radiation Fluxes

The space-integrated variability of the incident solar radiation is highly tied to the sun position, geographic latitude and terrain rugosity. Figure 8A,B show the effect of slope degree and aspect on the daily and instantaneous variability of basin-averaged incident solar radiation using an eight-direction classification system of slope aspect, during a typical, clear-sky day (e.g., 20 March 2016). To help interpret the topographic-induced variability of the cumulative diurnal radiation (Figure 8A), instantaneous radiation fluxes have been plotted against Local Standard Time (LST = GMT-7) per slope aspect, integrated over all the slope values (see Figure 8B). Due to the latitudinal location range of the basin (see Figure 1A) it is evident that the maximum energy flux is received by south-facing slopes (S, SW, SE), particularly those with a slope inclination between 42.77 and 43.24 , coincident with solar declination at the equinox for this basin, whose normal vector is parallel to the sun rays for most of this day. A unimodal pattern is observed in K i n for these aspects due to this effect. The maximum daily radiation received by south-facing slopes is K i n = 32.6 MJ/m 2 (Figure 8A) for which the largest direct heat flux, K i n = 1020 W/m 2 ), occurs between 12 and 1 PM LST (Figure 8B). Both SW and SE facing slopes reach similar incident values with slightly larger values given by the larger number of south-oriented elements in SW. A similar pattern is illustrated by both W and E facing slopes with maximum daily values of K i n = 24.2 MJ/m 2 , similar to those attained by upward-facing flat elements (F). NW and NE facing slopes reach their maximum daily radiation, K i n = 23.4 MJ/m 2 , for small slope inclinations with steady reductions for slope increments. Finally, the least illuminated surfaces are those facing north. Analogously to NW and NE, N facing slopes reach their largest radiation values for fairly flat inclinations and almost linearly reducing K i n until about 50 slopes. For slopes greater than 50 , K i n fluctuates around a rather low constant value of 4.5 MJ/m 2 . These evident differences in instantaneous and cumulative diurnal incident radiation that become greater as slope inclination increases are determinant in inducing contrasts in air and soil temperatures, prevalence of on-the-ground snow, evapotranspiration rates, surface soil moisture storage and vegetation phenology and type, commonly missing in coarse-scale land surface-atmosphere models.

4.3. Seasonal and Diurnal Variability of Incident Solar Radiation Components

For any particular day-hour (i.e., fixed sun position) the spatial variability of incident solar radiation is conditioned by the TIN element’s latitude, slope degree and aspect, topography shadow casting, sky and land view fraction, atmospheric transmissivities and land albedo. Contrastingly, for a particular TIN element (i.e., fixed geographic location) the incident daily radiation is controlled by the solar declination, atmospheric transmissivities and land surface albedo. Animations of the full day total incident radiation cycles are shown in Videos 1, 2, 3 and 4 for the spring equinox, summer solstice, fall equinox, and winter solstice respectively. At the daily scale, the diurnal range of energy flux at a TIN element, season-fixed, is controlled by the solar angle and the relative location of the solar vector and any remote shading. Figure 9 summarizes this spatio-temporal variability through day-accumulated maps of total incident radiation (K i n d a y in MJ/day) and diurnal distribution of total downwelling instantaneous solar radiation and its components at unshaded, flat-terrain TIN elements (in W/m 2 ). For purposes of comparison, clear-sky observations from the closest SURFRAD solar radiation station at Table Mountain Colorado (lon = 105.2368 W, lat = 40.12498 N, elevation = 1689 m) have been added to the diurnal plots organized by seasonal sequence. Use of a single, unshaded radiation station constitutes a limitation for model validation, particularly in narrow- valleys where remote shading and albedo-reflection is significant. The basin-wise incident radiation distribution maps (Figure 9A,C,E,G), on the left column, evidence the distinct amounts of solar incident energy per day-season ranging from an average of 34.5 MJ/m 2 in summer to 8.8 MJ/m 2 in winter. Despite the two equinoxes having the same solar declination, and thus energy input, basin albedo differs between the two ( α ¯ s p r i n g = 64% and α ¯ f a l l = 39%) due to the presence of on-the-ground snow in Spring that increases total incident daily radiation. A closer look to the diurnal cycles of computed and observed incident radiation in Figure 9B,D,F,H, reveal both similarities and differences with respect to observed fluxes from SURFRAD. Consistently across seasons, the peaks in K t o t a l s i m and K t o t a l o b s , k d i r s i m and k d i r o b s , and k d i f s i m and k d i f o b s show a slight delay, due to the time and longitudinal difference between the location of the SURFRAD station and the study catchment. Despite this fact, the simulated diurnal pattern illustrates the characteristic uni-modal distribution during the day with the corresponding seasonal variability. Additionally, the relative contribution of the diffuse radiation (k d i f ) is larger in summer, due to the increase in air temperatures that makes atmospheric transmissivities larger (see Appendix A). Finally, the contributions of K a l b and K b a c are larger during the Spring equinox, due to the higher influx of direct radiation and land albedo, but close to zero during the summer solstice when land albedo is heavily reduced. The maximum contribution of albedo reflected and atmospheric backscattered radiation occurs at noon during the Spring equinox when both contributions make up about 7% of the total incident radiation flux.

5. Discussion

The developed methodology and obtained results represent an improvement for the characterization of topographically-induced radiation differences in high-resolution, distributed models that use unstructured domain discretizations. The approach followed for the calculations of element slope and aspect, adds to the works of Evans [19], Hancock [20], Horn [21], Skidmore [22], Smith et al. [23,24], Zevenbergen and Thorne [25], Corripio [28] that calculated topographic properties based on regularly-spaced grid points or using a raster- (i.e., non vectorial) based algorithm. Furthermore, the methodology for estimating remote shading and sky view fraction adds to the state-of-the-art for its straightforward application in TIN meshes without need for refinement/derefinement procedures [30]. The use of horizon visibility-constraining distances, the consideration of difuse radiation, remote albedo, atmospheric backscattering and minimum algorithm inputs, constitutes an improvement respect to previous efforts by Marsh et al. [9], Clarke [31], Rinehart et al. [10] that partially consider direct and difuse components or that do not account for visible, sky-blocking features along all directions in irregularly-distributed node meshes. Despite its novelty, the algorithm and its application might have limitations regarding: (1) the accuracy of the original SRTM DEM; (2) the TIN sampling strategy and (3) the existene of a single solar radiation station in the neighboring area of study. The accuracy of the original DEM can introduce significant errors in the incident radiation flux calculation, particularly for the direct component, by missrepresentation of topographic features (e.g., slope and aspect) in rough topography and dense-canopy areas [49,50]. The TIN sampling strategy, to reduce the number of mesh nodes, affects the distribution of the solar radiation by providing accurate estimates in areas of high TIN density while leading to detrimental values in low-resolution areas [51,52,53]. We anticipate more accurate results in TINs with adaptive mesh densities in topograpically complex areas. Finally, the existence of a single solar radiation station limits our validation efforts at the catchment scale. The nationwide low-availability of radiation measurement sites (i.e., direct and diffuse components) is a constrain for our model effort. Finally, we recommend to test this algorithm in lager-size basins to evaluate its efficiency and accuracy relative to traditional grid-based methods. A further effort should focus on introducing the effect of vegetation and clouds in the calculations and quantification of the gains for hydrologic modeling in supercomputing, parallel applications.

6. Conclusions

The goal of this article was to present a TIN-based algorithm to compute slope, aspect and the clear-sky components of the incident solar radiation, accounting for direct, diffuse, atmospheric backscattered and albedo-reflected contributions that can be used for multiple applications, including hyper-resolution or large-scale land/atmospheric modeling. The method accounts for self-shading and cast shadows from remote, but prominent, topographic features and develops a pivot-based, distance-constrained method to compute sky view fraction with minimum information storage requirements. Using a TIN developed for this purpose and the definition of a planar normal vector, both slope degree and aspect maps were derived from the elevation TIN. Hourly weather information from a in-catchment weather station was used to compute optical airmass and atmospheric transmissivities and the NLDAS albedo product, at 12.5 km resolution, provided basin-averaged albedo values for each season. The algorithm was tested in a topographically complex basin in east-central Wyoming with observed data from a nearby SURFRAD solar radiation station. The spatio-temporal simulated outputs provide the following insights about the application of the model to this specific catchment:
  • In areas of high topographic complexity, like those on the east side of the Green River catchment, there are significant differences in instantaneous and daily cumulative incident solar radiation given by the high variability of slopes and aspects.
  • The method proposed here captures the seasonal variability of incident radiation and its diurnal cycle when compared with observed values.
  • The contribution of diffuse radiation is higher for the summer months due to the higher air temperatures and lower precipitable water.
  • The relative contribution of both atmospheric backscattered and landscape-reflected radiation is higher during the Spring season, when on-the-ground snow increases land surface albedo and thus more energy is received from the surrounding landscape and reflected back from the lower atmosphere.
Given the fast development and implementation of TIN-based models in hydrology and geomorphology [18,35,54,55], an improved scheme that helps resolve local differences in radiation, water and environmental dynamics is desirable. Scientifically and operationally, the methodology developed in this manuscript represents an improvement for the characterization of topographically-induced radiation differences in high-resolution, distributed models that use unstructured domain discretizations. Its direct application to TIN geometries will eliminate any further conversion from TIN to grid or the adaptation of fixed-orientation techniques to unstructuered data. This will significantly help with the understanding of fine-scale ecohydrologic dynamics otherwise difficult to represent with traditionally grid-based calculation methods. The code developed is free source and will be released by the authors upon request.

Acknowledgments

This research was funded by the Army Research Office award W911NF-18-1-0007, the NSF cooperative agreement EPS-1135483 and the NOAA/NWS Office of Water Predictions to the University of Wyoming. We are grateful with Robert Steinke from the Utah-Wyoming CI-Water project for providing the grid DEM for the small Green River catchment used in this article and OU MS student Jorge Celis for helping with figure-editing.

Author Contributions

H.A. Moreno, F.L. Ogden and L.V. Alvarez conceived and designed the model development and target experiments; H.A. Moreno and L.V. Alvarez developed code and figures; F.L. Ogden assisted with algorithm formulation and test design. H.A. Moreno led the paper writing.

Conflicts of Interest

The authors declare no conflict of interest.The founding sponsors had no role in the design of the study in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:
TINTriangular Irregular Network
DEMDigital Elevation Model

Appendix A. Solar Vector Calculation

This section elaborates on the procedure to estimate the solar vector from an orthogonal reference system with origin at the observer position. Please refer to Thomson [56] and Corripio [28] for a detailed description of the topocentric coordinate system. The variability of the solar declination can be obtained from different equations involving earth’s rotation and obliquity of the axis and other, less significant, earth’s movements, such as precession and nutation in the obliquity of the ecliptic. Bourges [57] provides a Fourier series approximation with a mean error of 0.008 and a maximum of 0.02 in terms of the day number D (see Equations (A1) and (A2)).
δ = 0.3723 + 23.2567 s i n D 0.758 c o s D + 0.1149 s i n 2 D + 0.3656 c o s 2 D 0.1712 s i n 3 D + 0.0201 c o s 3 D
where D is the day number:
D = ( 360 / 365.25 ) ( J 79.346 )
and J is the difference between the Julian day in consideration and the Julian day on January the first at noon of the given year, plus 1. At noon, the three coordinates of the solar vector in the east-west (s o x ), north-south (s o y ) and vertical (sphere-radial) directions (s o z ) are calculated using the element ϕ and δ , as shown in Equation (A3).
s o = ( 0 , s i n ( ϕ δ ) , c o s ( ϕ δ ) )
Having calculated the solar declination at noon, the estimation of s o , for any time (t), depends on the hour angle ( ω ) between the observer meridian and the solar meridian, by convention positive before noon and negative after noon. As the entire coordinate system rotates relative to solar noon, this movement can be decomposed in three rotations: one around the X-axis, to place the Z-axis parallel to the axis of rotation of the Earth; a second rotation around the Z-axis at an angle ω and a third rotation back around the X-axis to the observer position. To find the coordinates of the solar vector in the new, rotated reference system we multiply the original coordinates by three rotational matrices describing these movements, as shown in Equation (A3).
s = r x ( γ ) r z ( ω ) r x ( γ ) s o
with r being a rotation matrix around axis in subscript and angle in parenthesis, γ is the angle between Earth’s axis and the topocentric coordinate system Z-axis, and the hour angle ( ω ) is zero at noon and has the following value in radians at any time t (LST) given in hours and decimal fraction:
ω = π t 12 1
In matrix notation:
s = 1 0 0 0 c o s γ s i n γ 0 s i n γ c o s γ c o s ω s i n ω 0 s i n ω c o s ω 0 0 0 1 1 0 0 0 c o s ( γ ) s i n ( γ ) 0 s i n ( γ ) c o s ( γ ) S o x S o y S o z
Simplifying, the coordinates of the three components of the sun vector are:
s = s i n ω c o s δ s i n ϕ c o s ω c o s δ c o s ϕ s i n δ c o s ϕ c o s ω c o s δ + s i n ϕ s i n δ

Appendix B. Atmospheric Transmissivites

Atmospheric transmissivities ( τ ) can be calculated as indicated by Suckling and Hay [42] and Dingman [3]:
τ w a = ( 1 a w )
τ d a = ( 1 a d )
τ w s = ( 1 ρ w )
τ r s = ( 1 ρ r s )
τ d s = ( 1 ρ d s )
where a w and a d are the fractions of solar radiation absorbed by water vapor and aerosols (e.g., dust), respectively; and ρ w , ρ r s and ρ d are the fractions of solar radiation scattered by water vapour, air molecules (Rayleigh scattering) and dust and solid aerosols. These transmissivities are empirically related to optical airmass (M), humidity (RH) and air temperature (T) as:
τ w a = 1 0.077 ( M W ) 0.3
τ d a = 0.965 M
τ w s = 1 0.0225 M W
τ r s = 0.972 0.08262 M + 0.00933 M 2 0.00095 M 3 + 0.0000437 M 4
τ d s = 0.965 M
W is the precipitable water content in cm estimated as
W = 0.00493 R H T e x p 26.23 5416 T
where RH is the relative humidity (%) and T is the air temperature (K).
M is the optical air mass as suggested by Kasten and Young [58]. It has into account the curvature or Earth and is accurate to model the path thickness towards the horizon.
M = 1 c o s z + 0.50572 ( 96.07995 z ) 1.6364
At elevations different than the sea level, Equation (A19) needs to be corrected because of the smaller atmospheric thickness [3], which can be accounted for by calculating M(z).
M ( z ) = M e x p z 7000
where z is the elevation of each TIN element in meters.

References

  1. Beven, K.J. Rainfall-Runoff Modeling: The Primer, 2nd ed.; Wiley-Blackwell: Oxford, UK, 2012. [Google Scholar]
  2. Mirus, B.B.; Loague, K. How runoff begins (and ends): Characterizing hydrologic response at the catchment scale. Water Resour. Res. 2013, 49, 2987–3006. [Google Scholar] [CrossRef]
  3. Dingman, L. Physical Hydrology, 3rd ed.; Waveland Press, Inc.: Long Grove, IL, USA, 2015. [Google Scholar]
  4. Nassif, S.; Wilson, E. The influence of slope and rain intensityon runoff and infiltration. Hydrol. Sci. Bull. 1975, 20, 539–553. [Google Scholar] [CrossRef]
  5. Montgomery, D.; Dietrich, W.E. A physically-based model for the topographic control on shallow landsliding. Water Resour. Res. 1994, 30, 1153–1171. [Google Scholar] [CrossRef]
  6. Shen, H.; Julien, P. Erosion and sediment transport. In Handbook of Hydrology; Mc Graw Hill Professional: New York, NY, USA, 1993. [Google Scholar]
  7. Dozier, J.; Frew, J. Rapid calculation of terrain parameters for radiation modeling from digital elevation data. IEEE Trans. Geosci. Remote Sens. 1990, 28, 963–969. [Google Scholar] [CrossRef]
  8. Dubayah, R.; Rich, P.M. Topographic solar radiation models for GIS. Int. J. Geogr. Inf. Syst. 1995, 9, 405–419. [Google Scholar] [CrossRef]
  9. Marsh, C.; Pomeroy, J.W.; Spiteri, R. Implications of mountain shading on calculating energy for snowmelt using unstructured triangular meshes. Hydrol. Process. 2012, 26, 1767–1778. [Google Scholar] [CrossRef]
  10. Rinehart, A.J.; Vivoni, E.R.; Brooks, P.D. Effects of vegetation, albedo, and solar radiation sheltering on the distribution of snow in the Valles Caldera, New Mexico. Ecohydrology 2008, 1, 253–270. [Google Scholar] [CrossRef]
  11. Coop, J.D.; Givnish, T.J. Gradient analysis of reversed treelines and grasslands of the Valles Caldera, New Mexico. J. Veg. Sci. 2007, 18, 43–54. [Google Scholar] [CrossRef]
  12. Dozier, J. A clear-sky spectral solar radiation model for snow-covered mountainous terrain. Water Resour. Res. 1980, 16, 709–718. [Google Scholar] [CrossRef]
  13. Hood, E.; Williams, M.; Cline, D. Sublimation from a seasonal snowpack at a continental, mid-latitude alpine site. Hydrol. Process. 1999, 13, 1781–1797. [Google Scholar] [CrossRef]
  14. Garrote, L.; Bras, R.L. A distributed model for real-time flood forecasting using digital elevation models. J. Hydrol. 1995, 167, 279–306. [Google Scholar] [CrossRef]
  15. Aurenhammer, F.; Klein, R.; Lee, D.T. Voronoi Diagrams and Delaunay Triangulations; World Scientific Publishing Company: Singapore, 2013. [Google Scholar]
  16. Ivanov, V.Y.; Vivoni, E.R.; Bras, R.L.; Entekhabi, D. Catchment hydrologic response with a fully distributed triangulated irregular network model. Water Resour. Res. 2004, 40, W11102. [Google Scholar] [CrossRef]
  17. Slingsby, A. An Object-Orientated Approach to Hydrological Modelling using Triangular Irregular Networks. In Proceedings of the GISRUK03, London, UK, April 2003. [Google Scholar]
  18. Tucker, G.; Lancaster, S.; Gasparini, N.; Brass, R.; Rybarczyk, S. An object-oriented framework for distributed hydrologic and geomorphic modeling using triangulated irregular networks. Comput. Geosci. 2001, 27, 959–973. [Google Scholar] [CrossRef]
  19. Evans, I. An integrated system of terrain analysis and slope mapping. Z. Geomorphol. 1980, 36, 274–295. [Google Scholar]
  20. Hancock, G. The impact of different gridding methods on catchment geomorphology and soil erosion over long time scales using a landscape evolution model. Earth Surf. Proc. Landf. 2006, 31, 1035–1050. [Google Scholar] [CrossRef]
  21. Horn, B. Hillshading and the reflectance map. IEEE Proc. 1981, 69, 41–47. [Google Scholar] [CrossRef]
  22. Skidmore, A. A comparison of techniques for calculating gadient and aspect from a gridded digital elevation model. Int. J. Geogr. Inf. Sci. 1989, 323–334. [Google Scholar] [CrossRef]
  23. Smith, M.B.; Seo, D.J.; Koren, V.I.; Reed, S.M.; Zhang, Z.; Duan, Q.; Moreda, F.; Cong, S. The distributed model intercomparison project (DMIP): Motivation and experiment design. J. Hydrol. 2004, 298, 4–26. [Google Scholar] [CrossRef]
  24. Smith, M.B.; Koren, V.; Reed, S.; Zhang, Z.; Zhang, Y.; Moreda, F.; Cui, Z.; Mizukami, N.; Anderson, E.A.; Cosgrove, B.A. The distributed model intercomparison project—Phase 2: Motivation and design of the Oklahoma experiments. J. Hydrol. 2012, 418, 3–16. [Google Scholar] [CrossRef]
  25. Zevenbergen, L.; Thorne, C. Quantitative analysis of land surface topography. Earth Surf. Proc. Landf. 1987, 12, 47–56. [Google Scholar] [CrossRef]
  26. Chang, K.; Tsai, B. The effect of DEM resolution on slope and aspect mapping. Cartogr. Geogr. Inf. Syst. 1991, 18, 69–77. [Google Scholar] [CrossRef]
  27. Gao, J. Resolution and accuracy of terrain representation by grid DEMs at a micro-scale. Int. J. Geogr. Inf. Sci. 1997, 11, 199–212. [Google Scholar] [CrossRef]
  28. Corripio, J.G. Vectorial algebra algorithms for calculating terrain parameters from DEMs and solar radiation modelling in mountainous terrain. Int. J. Geogr. Inf. Sci. 2003, 17, 1–23. [Google Scholar] [CrossRef]
  29. Wood, E.F.; Roundy, J.K.; Troy, T.J.; van Beek, L.P.H.; Bierkens, M.F.P.; Blyth, E.; de Roo, A.; Döll, P.; Ek, M.; Famiglietti, J.; et al. Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water. Water Resour. Res. 2011, 47. [Google Scholar] [CrossRef]
  30. Montero, G.; Escobar, J.; Rodriguez, E.; Montenegro, R. Solar radiation and shadow modelling with adaptive triangular meshes. Sol. Energy 2009, 83, 998–1012. [Google Scholar] [CrossRef]
  31. Clarke, J. Energy Simulation in Building Design, 2nd ed.; Oxford University Press: Oxford, UK, 2001. [Google Scholar]
  32. Cano, D.; Monget, J.; Guillard, H.; Albuisson, M.; Regas, N.; Wald, L. A method for the determination of the global solar radiation from meteorological satellite data. Sol. Energy 1986, 37, 31–39. [Google Scholar] [CrossRef]
  33. Janjai, S. A method for estimating direct normal solar radiation from satellite data for a tropical environment. Sol. Energy 2010, 84, 1685–1695. [Google Scholar] [CrossRef]
  34. Zelenka, A.; Perez, R.; Seals, R.; Renne, D. Effective accuracy of the satellite-derived hourly irradiance. Theor. Appl. Climatol. 1999, 62, 199–207. [Google Scholar] [CrossRef]
  35. Ogden, F.; Lai, W.; Steinke, R. ADHydro-Quasy-3D high-performance hydrologic model. In Proceedings of the 3rd Joint Federal Interagency Conference (10th 46 Federal Interagency Hydrologic Modeling Conference), Reno, NV, USA, 19–23 April 2015. [Google Scholar]
  36. Mock, C. Climatic controls and spatial variations of precipitation in the western United States. J. Clim. 1996, 9, 1111–1125. [Google Scholar] [CrossRef]
  37. Vivoni, E.R.; Ivanov, V.Y.; Bras, R.L.; Entekhabi, D. Generation of triangulated irregular networks based on hydrological similarity. J. Hydrol. Eng. 2004, 9, 288–302. [Google Scholar] [CrossRef]
  38. Vivoni, E.R.; Mascaro, G.; Mniszewski, S.; Fasel, P.; Springer, E.P.; Ivanov, V.Y.; Bras, R.L. Real-world hydrologic assessment of a fully-distributed hydrological model in a parallel computing environment. J. Hydrol. 2011, 409, 483–496. [Google Scholar] [CrossRef]
  39. Iqbal, M. An Introduction to Solar Radiation; Academic Press: Toronto, ON, Canada, 1983. [Google Scholar]
  40. Lambert, J. Photometrie; W. Engelmann: Leipzig, Germany, 1760. [Google Scholar]
  41. Campbell, G.S.; Norman, J. An Introduction to Environmental Biophysics; Springer: New York, NY, USA, 1988. [Google Scholar]
  42. Suckling, P.; Hay, J. Modeling direct, diffuse and total solar radiation for cloudless days. Atmosphere 1976, 15, 194–207. [Google Scholar]
  43. Varley, M.J.; Beven, K.J.; Oliver, H.R. Modelling solar radiation in steeply sloping terrain. Int. J. Climatol. 1996, 16, 93–104. [Google Scholar] [CrossRef]
  44. Proy, C.; Tanre, D.; Deschamps, P. Evaluation of topographic effects in remotely sensed data. Remote Sens. Environ. 1989, 30, 21–32. [Google Scholar] [CrossRef]
  45. Wang, J.; Robinson, G.J.; White, K. A fast solution to local viewshed computation using grid-based digital elevation models. Photogramm. Eng. Remote Sens. 1996, 62, 1157–1164. [Google Scholar]
  46. Wang, J.; Robinson, G.J.; White, K. Generating viewsheds without using sightlines. Photogramm. Eng. Remote Sens. 2000, 66, 87–90. [Google Scholar]
  47. Nunez, M. The calculation of solar and net radiation in mountainuous terrain (Ridson, Tasmania). J. Biogeogr. 1980, 7, 173–186. [Google Scholar] [CrossRef]
  48. Xia, Y.; Mitchell, K.; Ek, M.; Sheffield, J.; Cosgrove, B.; Wood, E.; Luo, L.; Alonge, C.; Wei, H.; Meng, J.; et al. Continental-scale water and energy flux analysis and validation for the North American Land Data Assimilation System project phase 2 (NLDAS-2): 1. Intercomparison and application of model products: water and energy flux analysis. J. Geophys. Res. Atmos. 2012, 117, 1–27. [Google Scholar] [CrossRef]
  49. Rodriguez, E.; Morris, C.S.; Belz, J.E. A global assessment of the SRTM performance. Photogramm. Eng. Remote Sens. 2006, 72, 249–260. [Google Scholar] [CrossRef]
  50. Mukul, M.; Srivastava, V.; Mukul, M. Analysis of the accuracy of shuttle radar topography mission (SRTM) height models using international global navigation satellite system service (IGS) network. J. Earth Syst. Sci. 2015, 124, 1343–1357. [Google Scholar] [CrossRef]
  51. Shewchuk, J.R. Delaunay refinement algorithms for triangular mesh generation. Comput. Geom. 2002, 22, 21–74. [Google Scholar] [CrossRef]
  52. Vivoni, E.R.; Teles, V.; Ivanov, V.Y.; Bras, R.L.; Entekhabi, D. Embedding landscape processes into triangulated terrain models. Science 2005, 19, 429–457. [Google Scholar] [CrossRef]
  53. Vivoni, E.R.; Ivanov, V.Y.; Bras, R.L.; Entekhabi, D. On the effects of triangulated terrain resolution on distributed hydrologic model response. Hydrol. Process. 2005, 19, 2101–2122. [Google Scholar] [CrossRef]
  54. Francipane, A.; Ivanov, V.Y.; Noto, L.V.; Istanbulluoglu, E.; Arnone, E.; Bras, R.L. tRIBS-Erosion: A parsimonious physically-based model for studying catchment hydro-geomorphic response. Catena 2012, 92, 216–231. [Google Scholar] [CrossRef]
  55. Refice, A.; Giachetta, E.; Capolongo, D. SIGNUM: A Matlab, TIN-based landscape evolution model. Comput. Geosci. 2012, 45, 293–303. [Google Scholar] [CrossRef]
  56. Thomson, D.B. Introduction to Geodetic Astronomy; Department of Surveying Engineering, University of New Brunswick: Fredericton, NB, Canada, 1981. [Google Scholar]
  57. Bourges, B. Improvement in solar declination computation. Sol. Energy 1985, 35, 367–369. [Google Scholar] [CrossRef] [Green Version]
  58. Kasten, F.; Young, A.T. Revised optical air mass tables and approximation formula. Appl. Opt. 1989, 28, 4735–4738. [Google Scholar] [CrossRef] [PubMed]
Figure 1. General location of the uppper Green River basin, a major tributary to the Colorado River within the state of Wyoming in the United States of America. The figure also shows a DEM with elevations and stream channel networks.
Figure 1. General location of the uppper Green River basin, a major tributary to the Colorado River within the state of Wyoming in the United States of America. The figure also shows a DEM with elevations and stream channel networks.
Water 10 00398 g001
Figure 2. (A) Irregular-mesh model of study watershed; (B) Sample zoom region with TIN elements. Each triangle’s elevation is obtained from the mean of its three vertices’ values; (C) Normal vector to a triangular plane determined by the cross product of vectors u and v .
Figure 2. (A) Irregular-mesh model of study watershed; (B) Sample zoom region with TIN elements. Each triangle’s elevation is obtained from the mean of its three vertices’ values; (C) Normal vector to a triangular plane determined by the cross product of vectors u and v .
Water 10 00398 g002
Figure 3. Spatial distribution and central tendency statistics of the (A) n u x west-east and; (B) n u y south-north components of n u . Dark blue colors represent west and south facing triangles while dark red depicts east and north facing elements in A and B, respectively.
Figure 3. Spatial distribution and central tendency statistics of the (A) n u x west-east and; (B) n u y south-north components of n u . Dark blue colors represent west and south facing triangles while dark red depicts east and north facing elements in A and B, respectively.
Water 10 00398 g003
Figure 4. (A) Slope degree ( ζ ) and (B) slope aspect ( α ), for the Green River Basin. Slope aspects facing north-east are depicted in white while south-west facing slopes are represented black. Flat elements (zero-slope) are depicted in light grey.
Figure 4. (A) Slope degree ( ζ ) and (B) slope aspect ( α ), for the Green River Basin. Slope aspects facing north-east are depicted in white while south-west facing slopes are represented black. Flat elements (zero-slope) are depicted in light grey.
Water 10 00398 g004
Figure 5. Remote shading calculation in unstructured meshes. (A) Radial search for sky-blocking domain elements to each pivot (i.e., T 0 ) along n search slice directions (n = 16 in this example case). Dark and hollow circles represent included (higher) and discarded (lower) elements within each slice. The dotted lines represent the maximum horizon visibility, d v , and the grey-shaded areas extend to the sky-blocking element from T 0 . For a particular slice search direction (e.g., SSE) several neighboring elements above the horizon (e.g., z(T 0 ) < z(Ti) with i = 2 , 3 , 4 , 5 , 6 ) are included by testing that their horizontal azimuths fall between the slice limits (e.g., 67.5 and 90 ); (B) Δz/d Vs. d for the TIN elements within the SSE search slice. Note that T 3 is the sky-blocking element along this direction; (C) Zenith angle (f v s u ) calculation of the line of sight between the pivot T 0 and the sky-blocking vertex T 3 . Note that T 3 is within the observable horizon distance d v . The sun zenith angle ( ϕ δ ) is also calculated in terms of the distance, d, and the difference in elevation Δz. In this case since f v s u < ( ϕ δ ) then T 0 is shaded by T 3 .
Figure 5. Remote shading calculation in unstructured meshes. (A) Radial search for sky-blocking domain elements to each pivot (i.e., T 0 ) along n search slice directions (n = 16 in this example case). Dark and hollow circles represent included (higher) and discarded (lower) elements within each slice. The dotted lines represent the maximum horizon visibility, d v , and the grey-shaded areas extend to the sky-blocking element from T 0 . For a particular slice search direction (e.g., SSE) several neighboring elements above the horizon (e.g., z(T 0 ) < z(Ti) with i = 2 , 3 , 4 , 5 , 6 ) are included by testing that their horizontal azimuths fall between the slice limits (e.g., 67.5 and 90 ); (B) Δz/d Vs. d for the TIN elements within the SSE search slice. Note that T 3 is the sky-blocking element along this direction; (C) Zenith angle (f v s u ) calculation of the line of sight between the pivot T 0 and the sky-blocking vertex T 3 . Note that T 3 is within the observable horizon distance d v . The sun zenith angle ( ϕ δ ) is also calculated in terms of the distance, d, and the difference in elevation Δz. In this case since f v s u < ( ϕ δ ) then T 0 is shaded by T 3 .
Water 10 00398 g005
Figure 6. Sky view fraction fvs as calculated from the algorithm presented in Section 3.4.1.
Figure 6. Sky view fraction fvs as calculated from the algorithm presented in Section 3.4.1.
Water 10 00398 g006
Figure 7. Instantaneous incident solar radiation flux (K i n ) for the Upper Green River Basin during the spring equinox (20 March) in year 2016 at (A) 12 PM and (B) 5 PM LST, (GMT-7). Radiation flux maps account for direct, self and remote shading, diffuse, atmospheric backscattered, and albedo-reflected components.
Figure 7. Instantaneous incident solar radiation flux (K i n ) for the Upper Green River Basin during the spring equinox (20 March) in year 2016 at (A) 12 PM and (B) 5 PM LST, (GMT-7). Radiation flux maps account for direct, self and remote shading, diffuse, atmospheric backscattered, and albedo-reflected components.
Water 10 00398 g007
Figure 8. (A) Cumulative diurnal incident solar radiation (K i n ) in MJ/m 2 /day for different slope degrees and aspects and (B) hourly instantaneous incident radiation flux (K i n ) in w/m 2 for different slope aspects including flat TIN elements (F), during 20 March (i.e., spring equinox) of year 2016.
Figure 8. (A) Cumulative diurnal incident solar radiation (K i n ) in MJ/m 2 /day for different slope degrees and aspects and (B) hourly instantaneous incident radiation flux (K i n ) in w/m 2 for different slope aspects including flat TIN elements (F), during 20 March (i.e., spring equinox) of year 2016.
Water 10 00398 g008
Figure 9. Spatial distribution of simulated cumulative daily incident radiation (MJ/m 2 ) and diurnal cycle of simulated and observed incident radiation flux (w/m 2 ) and its components in the upper Green River Basin during a typical, clear-sky day of: (A,B) spring equinox, (C,D) summer solstice, (E,F) autumn equinox, and (G,H) winter solstice. The super indices sim (solid lines) and obs (dashed lines) stand for simulated and observed values. K t o t a l , K d i r , K d i f , K b a c and K a l b represent the total downwelling, direct, diffuse, atmospheric backscattered and albedo-reflected radiation components. K d b a = K d i f + K b a c + K a l b .
Figure 9. Spatial distribution of simulated cumulative daily incident radiation (MJ/m 2 ) and diurnal cycle of simulated and observed incident radiation flux (w/m 2 ) and its components in the upper Green River Basin during a typical, clear-sky day of: (A,B) spring equinox, (C,D) summer solstice, (E,F) autumn equinox, and (G,H) winter solstice. The super indices sim (solid lines) and obs (dashed lines) stand for simulated and observed values. K t o t a l , K d i r , K d i f , K b a c and K a l b represent the total downwelling, direct, diffuse, atmospheric backscattered and albedo-reflected radiation components. K d b a = K d i f + K b a c + K a l b .
Water 10 00398 g009

Share and Cite

MDPI and ACS Style

Moreno, H.A.; Ogden, F.L.; Alvarez, L.V. Unstructured-Mesh Terrain Analysis and Incident Solar Radiation for Continuous Hydrologic Modeling in Mountain Watersheds. Water 2018, 10, 398. https://doi.org/10.3390/w10040398

AMA Style

Moreno HA, Ogden FL, Alvarez LV. Unstructured-Mesh Terrain Analysis and Incident Solar Radiation for Continuous Hydrologic Modeling in Mountain Watersheds. Water. 2018; 10(4):398. https://doi.org/10.3390/w10040398

Chicago/Turabian Style

Moreno, Hernan A., Fred L. Ogden, and Laura V. Alvarez. 2018. "Unstructured-Mesh Terrain Analysis and Incident Solar Radiation for Continuous Hydrologic Modeling in Mountain Watersheds" Water 10, no. 4: 398. https://doi.org/10.3390/w10040398

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop