Fluvial Transport Model from Spatial Distribution Analysis of Libyan Desert Glass Mass on the Great Sand Sea (southwest Egypt): Clues to Primary Glass Distribution

Libyan Desert Glass (LDG) is a natural silica-rich melted rock found as pieces scattered over the sand and bedrock of the Western Desert of Egypt, northeast of the Gilf Kebir. In this work, a population mixture analysis serves to relate the present spatial distribution of LDG mass density with the Late Oligocene–Early Miocene fluvial dynamics in the Western Desert of Egypt. This was verified from a spatial distribution model that was predicted from the log-normal kriging method using the LDG–mass-dependent transformed variable, Y(x). Both low-and high-density normal populations (–9.2 < Y(x) < –3.5 and –3.8 < Y(x) < 2.1, respectively) were identified. The low-density population was the result of an ordinary fluvial LDG transport/deposition sequence that was active from the time of the melting process, and which lasted until the end of activity of the Gilf River. The surface distribution of the high-density population allowed us to restrict the source area of the melting process. We demonstrate the importance of this geostatistical study in unveiling the probable location of the point where the melting of surficial material occurred and the role of the Gilf River in the configuration of the observed strewn field.


Introduction
Libyan Desert Glass (LDG) is an entirely vitreous rock that occurs as fragments scattered within the sheet of coarse, poorly sorted, interdune (i.e., in "corridors" or "passages") sand or directly above the bedrock (i.e., in deflationary interdune areas).That is, LDG is located in flat areas between the north-south-trending longitudinal dunes of the Great Sand Sea in Libya and western Egypt at present.The LDG's age is about 29 million years [1][2][3]; it dates to the Oligocene epoch, when the region had a humid climate that produced intense precipitation.Runoff from the highlands (i.e., the Gilf Kebir Plateau) originated in the tributaries of the Gilf River [4,5] system, which flowed north along the actual Great Sand Sea, causing rapid weathering of strata and exposing layers of "Nubian Sandstone".
The first attempt to delineate the LDG-strewn field was made in 1932 [6].Spencer and Clayton underwent an expedition to reconsider the problem of LDG land distribution, and they proposed the first surface size value for the strewn field: approximately 80 km north-south and 25 km east-west [7].Spencer later indicated that the LDG extends over 137 km N-S and 56 km E-W [8].In 1984, Weeks et al. sampled 27 sites in some interdune passages of the Great Sand Sea, which were spaced 10 km apart from one another [9].Finally, Barakat et al. sampled LDG at 47 sites and noted that the majority of the LDG surface mass was found within two important accumulation zones [10].
The majority of studies that have studied the origin of LDG from a geochemical perspective [1][2][3][11][12][13][14][15][16][17][18][19][20][21][22][23] conclude or are based on the assumption that this glass is the result of a meteoritic [11,13,[15][16][17][18][19][20][21][22][23] or comet [24] impact process.On the other hand, the layered structure of the LDG pieces, as well as the lack of clear evidence of aerodynamic morphologies typical of splashform tektites and the existence of high temperature phases, are the principal empirical arguments that Boslough and Crawford [25], and Wasson [26], took into account for the validation of a model of formation of LDG based on a low altitude airburst simulation.Concomitantly, craters BP and Oasis have been unequivocally associated with impacts [27,28].The impact origin of another geoform named Kebira, which is located relatively close to LDG surface fragments and was identified by satellite observations [29], has not been confirmed [30].In the case of the crater-like structures included in the so-called Gilf Kebir Crater Field [31], prior studies demonstrated the hydrothermal origin of these structures [32]; this latter work marked the end of the debate launched in the work presented by Paillou et al. [33] in 2006, in which both impact and hydrothermal hypotheses were contrasted.
The main goal of this study is to demonstrate that the actual distribution pattern of LDG fragments on the Oligocene course of the Gilf River can be used to reconstruct the associated melting event.For this purpose, a distribution map of LDG surface mass is constructed by combining statistical univariate analysis with correlation and interpolation of georeferenced data.Surface mass data was compiled from the surveys of two of the previously mentioned authors [9,10] and the "Gilf Kebir 2007" expedition [34]; the latter was organized by one of the authors (Marius Ramirez) of the present study.Geomatic methods are intrinsically associated with geostatistical work; they will be used to determine the paleodrainage network at the regional and local levels so as to establish the role of paleohydrography in the current distribution of LDG.

Study Area
Three geographically partitioned zones represent the interrelated scenarios investigated in this study: the Gilf River paleo-course (22°-32° N and 23°-29° E), the LDG strewn field (24°-26° N and 25°-26° E; insert in Figure 1), and its proximal geological and hydrographic context (24°-26° E and 24°-26° N).The first includes the last two and the study area at the regional level; it extends from the highlands of the Gilf Kebir plateau to the Mediterranean Sea along the Egypt-Libya border and represents the limits of our general paleohydrological study.Local-level study of Gilf River tributaries and the drainage pattern of the Kebira structure (suggested here as a product of fluvial erosion) will imply further partitions of the regional area.The proposed margins of the strewn field zone ensure the inclusion of all sampling points and appropriate resolution of the distribution map. Figure 1 shows the margins of every zone.The LDG strewn field's extent was maximized to display the relevant information supplied by the different surveys of the region.An area that encloses the mentioned strewn field and its western extension is also showed.This zone represents a second level of spatial study where relevant elements-for example, the connection between the Kufrah and Gilf River basins, or the accepted BP and Oasis impact craters-are included.This study partition is named "Extended Strewn Field" in Figure 1 and the goal with this zone is to establish the role of some wadis from the Gilf Kebir Plateau, as well as of relevant geoforms, in the formulation of a hypothesis about the origin of the LDG.
In the western part of the LDG-strewn field and the southwestern part of the Great Sand Sea, the Gilf Kebir Plateau arises 300 m above the desert floor and represents a topographical barrier between the Oligocene Gilf River and the Miocene Kufrah basins.It is an 8000-km 2 flat plain with numerous wadis; some of them, in the eastern flank of the plateau, conform to the ancient tributary system of the Gilf River.The northern part of this plateau is known as the Abu Ras Plateau, which is close to the LDG-strewn field.The Great Sand Sea is mostly composed of longitudinal sand dunes that extend for tens of kilometers; as high as 100 m, they have a quasi-north-south parallel arrangement.The northern depressions (Siwa and Qattara) are relics of fan deltas from the Oligocene (Figure 2), from the beginning of the Oligocene until the early-Miocene Epoch; the gradual uplift of southwest Egypt produced the regression of Tethys Sea northward, setting the shore to the Siwa or Qattara latitudes.The depressions were formed by the running water of the Gilf River that reached the shore, carving the depressions by dissolution of carbonates that were underlying the more resistant Miocene limestone dipping-northward-layers.This latter hypothesis on the formation of the northern depressions of the Western Desert of Egypt was supported by the discovery and study of paleorivers detected by radar, and it was formulated in conjunction with the fact that a humid climate during long time periods favored the recharge of the Gilf River from the highlands of the Gilf Kebir Plateau [4,35].The study area and the different reported Libyan Desert Glass (LDG)-strewn fields.Brown ellipse: the smallest distribution zone [6]; red line: the largest distribution zone [9]; orange line: the common western limit of LDG deposition [8].The bedrock materials underlying the aeolian sands and dunes are included in the "Nubian Sandstone".This name refers to the Nubian (Late Jurassic-Lower Cretaceous) depositional cycle; it is not a formal stratigraphic unit [36].In the study area, this continental sequence consists of the Gilf Kebir Formation [37], which is composed mainly of red quartzitic sandstones and conglomerates; on the higher parts of the Gilf Kebir Plateau, it is overlaid unconformably by kaolinitic sandstones of the Sabaya Formation, the latter of which are of Upper Cretaceous age [38].The Maghrabi Formation conformably overlies the Sabaya Formation and is composed of fine-grained sandstone that covers some parts of the Gilf Kebir Plateau (stratigraphic column in [38] page 311).The formations of "Nubian Sandstone" rest unconformably on sandstones and glaciofluvial diamictites of Paleozoic age (i.e., the Wadi Malik and North Wadi Malik Formations, respectively) in the Abu Ras Plateau.All these formations are possible candidates for having become LDG precursor target materials.Further additions to this group are the gneisses of Precambrian age, which have isotopic composition compatible with LDG [3], even though they are located in Gebel Oweinat, the "triple-point" of the borders of Egypt, Libya, and Sudan, which is approximately 250 km south of the LDG-strewn field.

Sampling and Pattern of Univariate Data
The stationary variable [39] Z(x) studied here is the mass of LDG distributed on the Great Sand Sea's surface.We assume that the total initial mass of LDG is homogeneous at the compositional level and that the degree of fragmentation is a negligible index of the surface distribution.Geostatistical study of the regionalized variable has remarkable limitations at the sampling stage, as accessing the sampled places is difficult because of hyperarid weather conditions and the presence of dunes.This difficulty was attenuated through the use of data collected during four different expeditions: "Gilf Kebir 2007" [34], Weeks et al. [9]-these authors reported two expeditions-and Barakat et al. [10].Thus, we conducted random sampling whereby 76 sampling points were provided by the only four known expeditions that focused on the mass distribution of LDG.The sampling points' locations were considered as the central ones within each area where LDG fragments were collected; thus, the variable was standardized dividing by the sampling area (Z(x)/Area (A) (g/m 2 )).The source of every sampling point is explained in Table 1 of the Supplementary Material.The sampling design was not based on a grid or specified terrain divisions; thus, it was not straightforward to derive an original extension of the strewn field.
The complete spatial randomness (CSR) [40] of the mapped point pattern is evaluated by contrasting the nearest neighbor index (NNI or R) [41] at small scales of the point process, with the graphical assessment of the K and L functions [42] for small and large distances between pairs of points.
Using the paleontological statistics (PAST) software package [43], the nearest neighbor analysis follows [41].The area is estimated by the smallest enclosing rectangle (Figure 3, 4.17 × 10 3 km 2 ).The null hypothesis is a random Poisson process (H0: events exhibit complete spatial randomness), giving a modified exponential nearest neighbor distribution with mean: where A is the area and n the number of sampling points.The probability that the distribution is Poisson-type depends on the R value:

Exploratory Statistics
Descriptive statistics were performed to establish measures of central tendency and dispersion.The non-normal (i.e., logarithmic) distribution of the values is evident in Figure 5a, in which the histogram has a large right tail (and some outliers among the greatest values).This distribution can be straightforwardly described by estimations of robust statistical parameters, such as the median (0.27 g/m 2 ) and the interquartile range (8.48 g/m 2 ) for central tendency and dispersion, respectively.The value of the asymmetry coefficient, or skewness (3.3595), is markedly positive-that is characteristic of a logarithmic distribution.In a log-normal distribution (Figure 5a), high values increase the variance, which increases the difficulty of interpretation, calculation of the variogram, and the subsequent kriging interpolation.For this reason, we transformed the standardized variable Z(x)/A (g/m 2 ) to Y(x) = ln[(Z(x) + 1)/A] (g/m 2 ), where the constant 1 was introduced to avoid an indeterminate mass value (i.e., ln0).This transformed variable's histogram suggested the existence of two modes (Figure 5b).Subsequent analysis of the mixture of populations was performed using the expectation-maximization (EM) algorithm [44] implemented in the PAST software package [43].
This analysis produced probabilities of ≈71% and ≈29% of finding a sampling point within the high-density (-3.8 < Y(x) < 2.1) and low-density (-9.2 < Y(x) < -3.5) populations, respectively.Table 1 shows the statistical features of both populations: means, standard deviations, and the probability p of belonging to high or low-density populations.
In order to determine the original population (and both of its sub-populations) of the mass per square meter, as well as its logarithmic transformation, two normality tests were performed by the PAST software: the Shapiro-Wilk test [45] and the Jarque-Bera test [46].
For the two tests the null hypothesis (H0) is as follows: the sample was taken from a population with normal distribution.If the given p (normal) is less than 0.05, then normal distribution can be rejected.The results are summarized in Table 2.For the original distribution, as well as for its logarithmic transformation, the value of p (normal) is less than 0.05; however, when the subpopulations are considered, the values of p (normal) are above 0.05 and conditions of normality are accepted.In the case of the Shapiro-Wilk test, normality was accepted with a Wcrit (critical value of the Shapiro-Wilk test) of 0.90 and a probability-p (normal)-greater than 0.05.The Jarque-Bera method was in agreement with the results from the Shapiro-Wilk procedure.

Kriging Methods
The semi-variogram is a quantitative descriptive statistic used to characterize the spatial continuity of a regionalized variable and is represented by the function: where, γ(h) = value of semi-variogram; Y(xi) = value of log-transformed variable at spatial location (xi); and N(h) = number of pairs with separation distance h.
In order to compute all the pairs of sampled points, we use a polar grid whose cells are defined by a lag distance and an angle.Only pairs falling into a grid cell were computed for calculation of the variogram.The value of h is calculated from the differences between the angles and the lag values of the respective points in the pair.A lag width of 5000 m is selected because it approximates the average distance between neighboring sample points.Using Equation (3), we plot the omnidirectional semi-variogram of the variable ln[(Z(x) + 1)/A] together with directional semi-variograms oriented at 15°, 30°, 45°, 60°, 75°, 90°, 105°, 120°, 135°, 150°, and 165°.We chose the Gaussian mathematical model, since it agreed best with the data after fitting its parameters (sill, nugget effect, and range); then, we modeled the data's anisotropy with the directional semi-variograms.
We used the information from the semi-variogram analysis in the log-normal kriging to achieve optimal weighting functions, since ordinary kriging was applied to the logarithmic transformation of the data.The log-normal predictor kriging is [47]: Because the interpolation was realized in the variable Y(x) = ln[(Z(x) + 1)/A], the data must be back-transformed to the original measurement scale (g/m 2 ) by an adaptation of Yamamoto's method [48]: where, ( ) is the smoothing effect described by Yamamoto [49].
The reliability of any interpolation kriging should be checked by examination of its errors.A common method of error examination is residual estimation, or calculation of the coincidence degree between the interpolation and the original data.The residual is the difference between the sampled value and the value of Z interpolated in the same location inside the surface grid.The formula used to determine the residuals is: = predicted value at the same location.
In addition to residuals we calculate the mean error (ME) and the root mean square error (RMSE) to measure the bias and precision of kriging: ( ) ( )

Digital Tools on Image Processing
The Western Desert paleodrainage network and the drainage pattern [50] of the Kebira structure were constructed according to the D8 algorithm [51] using a digital elevation model of Shuttle Radar Topography Mission (SRTM) images [52].The paleodrainage network was constructed in a Geographic Information System (GIS).Moreover, extensions "basin 1" [53] and "Hydro" [54], both based on algorithm D8, were also used.The process consists of the following three steps: 1. Preparation of the Digital Elevation Model (DEM).
a. Elimination of the sinks.In order to obtain a "corrected DEM" for subsequent analysis, this function identifies depressed zones and fills them.b.Flux direction.The algorithm D8 traces the flow from every pixel of the DEM to one of the surrounding pixels.The result is the generation of flux lines that appear as lines perpendicular to the elevation contours.2. Determination of a raster with flow accumulation within each downslope cell.It is a calculation of the number of cells that flows into a particular downslope cell, according to flux direction [54].3. The paleodrainage network is a shape that results from the combination of the flux direction and accumulation layers by means of extension basin 1 [53].
Both paleodrainage maps served as base geographical guides to form the overlapped distribution LDG mass map and to evaluate the subsequent hypotheses about the possible points of LDG formation, as well as the directions and effects of fluvial transport.

Geostatistical Mapping
The omnidirectional experimental variogram was adjusted to a theoretical Gaussian model (Figure 6).Subsequently, directional variograms were adjusted following the same model, showing different vertical scales in their respective directions.This behavior indicates that the model has a geometric anisotropy (Figure 7).In the anisotropy direction of 30° the apparent length and the scale of the experimental variogram is maximized.This latter variogram is consequently selected for adjusting to the Gaussian plus nugget model (Figure 8).The adjusted parameters of the omnidirectional and the selected directional final model are summarized in Table 3.
The nugget effect is the discontinuity of the semi-variogram at the origin; both the omnidirectional and final directional semi-variograms present positive nuggets.In this study, we consider the nugget effect as a sum of variabilities at small distances, as well as the effects of the errors in sampling and position estimation.The range corresponds to the "influence zone" of the studied variable: at distances beyond the range value, the pairs of sampled points are no longer correlated.The final directional semi-variogram presents more spatial correlation among the sampled points (range = 20 km) than the omnidirectional semi-variogram (range = 17 km); this indicates the variation by direction of the LDG distribution in the Great Sand Sea.
The sill is the limit value at which the semi-variogram becomes more or less stable; it can be conceptualized as the variance of a regionalized variable and expresses the degree of the relationship between the sampled points on the studied surface.In the discussed semi-variogram models, the fact that sill values vary by direction (Figure 6 and Table 3) shows geometric anisotropy.We model the anisotropy at 30° because the vertical scale is maximized in that direction, so the semi-variance and the degree of relationship between the sampled points are also maximized.The direction of the final semi-variogram model is 120° (i.e., orthogonal to the direction of the maximum vertical scale).
The final semi-variogram model (Figure 8) was used to construct the contour map from the interpolation kriging, but we transformed the variable ln[(Z(x) + 1)/A] back to the original measurement scale (g/m 2 ).The corresponding map is illustrated in Figure 9, beside the distribution of the variable ln[(Z(x) + 1)/A] and the maps of the residuals (i.e., the true error of log-normal kriging).
We estimated the ranges of certain, probable, and possible anomalies by combining the means and standard deviations of the data (Table 1).For the low-density population, the lowest observed value (0.00010) fell within the range of certain anomalies (between 0.00004 and 0.00017), near the lower limit of the range of probable anomalies.This means that an almost-void value of LDG mass per square meter is a certain anomaly.
Within the high-density population, there was one value in the range of certain anomalies (between 8.214 and 34.535).This value indicates a peak of LDG accumulation on the surface.The limits of anomalies serve to illustrate the contour ranges on the distribution map.
The residuals map measures the true error of kriging.True errors less than and greater than zero indicate overestimated and underestimated interpolated values, respectively.The southern region may have more overestimated values, corresponding with the less-sampled area; in contrast, the area with zero error coincides with the accumulation of more sampling points.In addition to residuals, we calculated the ME and the RMSE.The ME is equal to −0.56 g/m 2 , showing that the log-normal kriging is an unbiased estimator.The average magnitude of the forecast error, measured by the RMSE indices, equals 4.65 g/m 2 ; it is worth mentioning that the RMSE is influenced more strongly by large errors than small ones.

Understanding the Surface Distribution Model of LDG
The Kufrah basin fluvial system is observed in the western part of the paleodrainage map (Figure 10), with its tributary dendritic system on the western edge of the northern Gilf Kebir.This system was active from the mid-Miocene to the Holocene [55].The present study, in comparison with earlier reports, demonstrates that the outfall of this Miocene river was probably at latitude farther north and west of Siwa.Similarly, this assessment validates the Siwa's depression as a part of the Oligocenic shore where the Gilf River would outflow [35].The map clearly illustrates the role of the Gilf Kebir elevation as the recharge zone of all the fluvial systems that emptied on Tethys at Qattara and Siwa latitudes and (specifically for the Kufrah basin) in northern Cyrenaica (an emerged zone from the Eocene).One limitation of determining flux direction from SRTM data is the erroneous identification of passages between dunes as subparallel drainage courses.This dune system is interpreted as aeolian deposits that developed over the fluvial course in dry epochs from the cessation of the streamflow (i.e., similar to the case of the Namib Desert in Namibia).The current aeolian dynamic differs slightly from the flux directions of Oligocene fluvial system, though a common and predominant south-north direction is evident (i.e., from the south to the north, and vice versa, for the fluvial and aeolian dynamics, respectively).Thus, in the scenario of an LDG deposition process from fluvial transport, the southernmost LDG fragments were likely deposited first.Figure 10b is a maximized selected view of the paleodrainage network of Figure 10a.It shows that one of the paleochannels with a predominant north-south-trending flow (the dendritic Wadi Abd el Malik; course represented by light blue line in Figure 12) crosses the Kebira sub-basin on its central and eastern part.Although this wadi has been understood as a tributary of Kufrah's Miocene basin [56], or even an independent fluvial system [29], it presents a certain probability of connection with the Gilf River fluvial system.The disposition of Wadi Abd El Malik (mainly for its medium course) is closed to the sources of the Gilf River tributaries; thus, at some moment, both systems could have been interconnected at this transition zone by a surface process of headward erosion.Moreover, the SRTM data in this zone indicate that both Wadi Abd el Malik and the sources of the tributaries of the eastern basin (Gilf River system) are situated at approximately the same elevation.Figure 10b shows an elliptic area illustrating this transitional connection zone between both systems.This scenario reflects how the Oligocene basin, east of the Gilf Kebir Plateau, contained tributaries from higher parts, including those within the Wadi Abd El Malik system (Figure 10b); subsequently, during the Miocene, some of them were assimilated into the Kufrah system of fluvial tributaries, as indicated by our paleodrainage model (Figure 10a).This is an example of west-east-trending dendritic-type highland connection between two principal fluvial systems; this area is located at the south of the southern flank of the Wadi Qubba, which is a wide, buried paleochannel with a slightly positive lenticular topography.This paleochannel is not inferred in the drainage map because of its moderate slope toward the northeast, which suggests that it was a zone of braided streams favoring high-volume sediment deposition, not a clear example of an erosive agent [56].Its southwest-northeast direction also indicates a Miocene stream connecting both principal basins in the region, but in this case it would have acted as a low-flow tributary of the Gilf River.
A paleodrainage study of the previously reported circular geoform named Kebira was conducted in order to obtain further criteria to uphold its possible impact origin.Although previous studies, and our observations while surveying the region, demonstrated a lack of petrological and mineralogical impact traces, Figure 11 shows how the Kebira structure contains all the types of drainage (i.e., wadis) that Mihályi et al. associated with impact structures [51].We assume that the hydrological pattern may help to identify ancient impact craters, but this is not a diagnostic criterion for identifying impact structures.Thus, the Kebira structure will remain as a simple geographical reference point during the discussion of the paleodrainage systems' behavior and the relative positioning of confirmed impact structures in the zone.
From the distribution map, we propose the differentiation of four distribution zones (Figure 12): Zone 1: among the two high-density LDG population zones identified, this is the northernmost.It includes the majority of the sampling points of normal population 2 and the largest pieces of LDG.Its boundaries roughly delineate an oval shape encompassing 25.2°-25.6°N and 25.5°-25.7°E. All the samples in this area were located over the interdune corridors, where some sandstone and siltstone outcrops of Upper Cretaceous formations are recognizable.The error in the kriging prediction (i.e., the residuals between sampled and predicted values) is very small in this zone, because the great number of sampling points in a relatively small area yielded a variogram model that indicated high correlation at small distances.Although this zone includes peak density values, its aggregation-trending pattern and the significant presence of data points with medium values determined a specific predicted map dominated by points that were not among the highest density values.Zone 2: this zone denotes the spatial correlation of both of the high-density LDG population zones at their northern and southern boundaries.At its western and eastern limits, a lack of LDG is predicted.Indeed, this suggests that there is confinement of the transported LDG along the Gilf River's path.This paradigmatic zone represents effective LDG transport by the Gilf River in the actual strewn field (purple arrow in Figure 12).Zone 3: this is the southernmost high-density LDG population.All the sampled points belong to population 2 (Table 1), but their scarcity produces large residuals during the prediction process via the log-normal kriging method (Figure 8).The location of this zone at the southernmost part of the strewn field, in addition to the existence of zone 2 (which functions as a transitional area between zone 3 and zone 1 to the north), lead us to consider it as the original depositional zone.Thus, the LDG fragments in this zone have suffered no (or little) effective fluvial transport.
Zone 4: the last partition is the northernmost zone; it contains an important predicted accumulation, though the errors are great at its western boundary.The mass isolines at zone 4 trend east-west along its long direction; this is its main difference from the other areas, which have a long direction with a dominant north-south component.This difference suggests an accumulation with a different origin from those of zones 1 and 3; furthermore, zone 4 has no southern transition zone like zone 2. Zone 4 could be associated with the sediment band of a meander bend of the Gilf River or a confluence of tributaries or streams.
Zone 1 is located adjacent to the southwestern endings of the Great Sand Sea (formally, a section where there are no longitudinal dunes), within a probable alluvial fan related to wadi Qubba (Figure 13).This is an east-west paleochannel that connects the Kufrah and Gilf River basins (i.e., it plays an equivalent role to the Wadi Abd El Malik in the uppermost parts of the Gilf Kebir Plateau).LDG fragments are seen in zone 1 because of the two principal dynamics that form the present desert landscape of the eastern Great Sand Sea: the erosive work of the Gilf River and the aeolian system that subsequently shaped the fluvial sands into dunes.In other words, the evident high concentration of LDG fragments in zone 1 is likely a consequence of the erosive work of the Gilf River over the LDG-bearing materials deposited in channel mouths (i.e., the alluvial fan) from the Wadi Qubba system.While LDG (which originally would have occurred as clasts into a sandy matrix, forming a sandstone or conglomerate) remained in situ after the erosion process, the sandy part of the sediment contributed to dune formation and sand sheets during dry periods.
Wadi Qubba is practically fully buried by Miocene and Quaternary sediments, which prevent or inhibit the outcropping of the Nubia Group or Carboniferous or Devonian materials, and subsequently, the LDG.The only area where we consider infilling materials to have been eroded is zone 1, where LDG fragments are scattered along the interdune surfaces and represent a distal part of an alluvial fan (or a flood plain).Thus, the Wadi Qubba probably acted as a tributary of the Gilf River.A line connecting the central location where the kriging method predicted the absence of LDG fragments (in the western part of zone 2) with the center of zone 1's peak accumulation follows the Qubba's path to the northeast (i.e., it is coincident with the Kufra-Abu Mingar camel route, mapped in 1939 [8]).In addition, this north-trending line marks the maximum gradient in LDG surface mass density (red arrow in Figure 12).
The LDG distribution is represented in Figure 13, where an accumulation is differentiated in the southern section.This feature (Zone 3) is marked as the original strewn field and its origin becomes compatible with a low-altitude airburst, i.e., the ground point contact of the possible fireball fall into Zone 3 of our distribution model.This scenario implies a subsequent regular transport-deposition of glass all along the south-north path of the Gilf River, from Zone 3 (the source of melted material) to Zone 4 (Figure 12).This hypothesis is reinforced if the discovery reported in [57] is considered; these authors collected samples of melt rock without high-pressure phases among sandstone layers about 160 km south-southeast of the LDG area (in the words of the authors).It is interpreted that this location refers to the southern part of Zone 3 because it is usual to consider the "LDG area" as our Zone 1, approximately.Its chemical composition suggests a melting temperature above 1600 °C.
The central accumulation (Zone 1 in Figure 12) is located northeast of the Saad Plateau and is traversed by the tributaries arising from the area of Wadi Qubba (Figure 13).The fragments in the northern area were transported from the maximum accumulation zone.This scenario explains why the area is considered to be formed by the mixture of both populations identified in the geostatistical study.
The estimation of the total surficial LDG mass was performed directly from the contour map of g/m 2 (Figure 13).In order to avoid overestimation of the mass of population 2, a "mixture" subpopulation associated to the more extended area (i.e., it is attributed to contours with intermediate values in Figure 13) is established (Table 4).Calculation of the total LDG mass resulted in a very low value compared with the simulation made by Artemieva et al.In the case of the Ries-Moldavite strewn field [58], they obtained the quantity of 1.6 × 10 13 g of tektites deposited in the simulation; on the other hand, Artemieva estimates a total mass of ejecta of ~ 4 × 10 15 g for the Australasian tektite strewn field [59], but for the LDG-strewn field Weeks et al. [9] calculated a mass of 1.4 × 10 9 g, while Barakat et al. [10] estimate a LDG mass equal to 2.67 × 10 8 g.Concerning these latter works, the result in Table 4 is very similar to the former work and almost one order of magnitude greater than the second.The occurrence of a layer of buried melted material is not ruled out in the Zone 3. In this sense, a radiative melting of target material from a low-altitude airburst is also compatible with the shape of the southernmost zone in the strewn field.

Conclusions
A general description of the transport process of Libyan Desert Glass fragments has been established through a geostatistical approach that comprises the structural spatial analysis of the sampled points and the kriging method.A great oval zone of important fragment accumulation is associated with the original source of LDG by a melting process (low altitude airburst or impact point without a visible crater) in the southern part of the strewn field.A second accumulation zone is separated from the former one by a transitional area with a small LDG mass surficial density that is associated to the Gilf River pathway with a dominant transport process, without clear evidence of massive accumulation.This second accumulation zone probably originated by a differential erosive action in this part of the Gilf River course or, alternatively, through the presence of a fluvial pool or lagoons in this section of the Oligocene/Miocene river.Similarly, we do not rule out the existence of a contribution of LDG mass from western tributaries in order to justify the great mass of accumulation in the mid-course of the Gilf River.In general, trends in geostatistical analysis marked the pathway of a Oligocene/Miocene river.This study demonstrated the usefulness of statistical analysis to elucidate the fluvial transport of mass from little evidence of differentiated melted material present on the surface of the Western Desert of Egypt.

Figure 1 .
Figure1.The study area and the different reported Libyan Desert Glass (LDG)-strewn fields.Brown ellipse: the smallest distribution zone[6]; red line: the largest distribution zone[9]; orange line: the common western limit of LDG deposition[8].

Figure 2 .
Figure 2. Principal geomorphological features in the Western Desert, Egypt.

2 )
where d is the observed mean distance between nearest neighbors.Clustered points produce R < 1, Poisson patterns give R ≈ 1, and overdispersed (regularity) points give R > 1. K(d), L(d), and L(d)-d functions represented in Figure4indicate a clustered pattern of points because the observed functions are above the high envelope (i.e., above the line of the 0.95 quantile) for distances up to 18 km.Thus, a clustered pattern of points is observed for short and large distances.

Figure 3 .
Figure 3.The smallest rectangle enclosing the LDG sampling points has been drawn over the geographic map.The rectangular area (4.17 × 10 3 km 2 ) was employed in the nearest neighbor analysis (R = 0.87), giving a clustered LDG point pattern.

Figure 4 .
Figure 4.The function for the observed LDG distribution appears in black; the three graphics correspond to functions K(d), L(d), and L(d)-d, which are above the 95% confidence interval (red functions), indicating a tendency towards a clustered distribution.

Figure 5 .
Figure 5. (a) Histogram of the variable g/m 2 of the LDG sampled points; (b) Histogram of the logarithmically transformed variable.
measured value of Z at location xi; and ( )

Figure 6 .
Figure 6.Omni-directional semi-variogram shows adjustment of the theoretical model to the experimental curve.Blue line is the mathematical semi-variogram model; the estimated semi-variogram is represented by the black line.

Figure 10 .
Figure 10.Paleodrainage maps: (a) Western Desert of Egypt; (b) inset of the red square in a, where the red ellipse indicates a probable connection zone between the Kufrah and Gilf River basins.

Figure 11 .
Figure 11.Drainage pattern of the Kebira structure.Structural units of the Kebira were proposed in [23].

Figure 12 .
Figure 12.Map of the LDG distribution model.Red arrow: line of maximum gradient in LDG surface mass density.Purple arrow: line of LDG transport.

Figure 13 .
Figure 13.Relation between paleohydrography and contour map of LDG distribution over the Landsat-7 ETM+ reference map.The red dotted line delimits the most likely area associated to the melting event.

Table 1 .
Data of the mixture populations (transformed variable).