Characterization of Black Sand Mining Activities and Their Environmental Impacts in the Philippines Using Remote Sensing

Magnetite is a type of iron ore and a valuable commodity that occurs naturally in black sand beaches in the Philippines. However, black sand mining often takes place illegally and increases the likelihood and magnitude of geohazards, such as land subsidence, which augments the exposure of local communities to sea level rise and to typhoon-related threats. Detection of black sand mining activities traditionally relies on word of mouth, while measurement of their environmental effects requires on-the-ground geological surveys, which are precise, but costly and limited in scope. Here we show that systematic analysis of remote sensing data provides an objective, reliable, safe, and cost-effective way to monitor black sand mining activities and their impacts. First, we show that optical satellite data can be used to identify legal and illegal mining sites and characterize the direct effect of mining on the landscape. Second, we demonstrate that Interferometric Synthetic Aperture Radar (InSAR) can be used to evaluate the environmental impacts of black sand mining despite the small spatial extent of the activities. We detected a total of twenty black sand mining sites on Luzon Island and InSAR ALOS data reveal that out of the thirteen sites with coherence, nine experienced land subsidence at rates ranging from 1.5 to 5.7 cm/year during 2007–2011. The mean ground velocity map also highlights that the spatial extent of the subsiding areas is 10 to 100 times larger than the mining sites, likely associated with groundwater use or sediment redistribution. As a result of this subsidence, several coastal areas will be lowered to sea level elevation in a few decades and exposed to permanent flooding. This work demonstrates that remote sensing data are critical in monitoring the development of such activities and their environmental and societal impacts.


Introduction
There are two types of iron ore most commonly used in the production of steel worldwide: hematite (Fe 2 O 3 , 69.9% Fe) and magnetite (Fe 3 O 4 , 72.4% Fe).Hematite deposits have been easily accessible and largely exploited until the last decade.Depletion of these deposits led to an increased demand for magnetite as the need for high-quality iron products continues [1,2].In the Philippines, magnetite occurs naturally in black sand, which results from the weathering and erosion of metamorphic and igneous rocks.Black sand accumulates in streams and drainage systems and is carried and deposited onto beaches.As a result, much of the Philippine coast is composed of black sand beaches.In response to the demand for magnetite, black sand mining and processing activities have significantly increased in recent years and the extracted magnetite is largely exported to China's steel mills [2,3].
Black sand mining occurs both legally and illegally but monitoring presents significant technical challenges and social costs.Reports from local civil service organizations [4] and the recent indictment of a provincial governor on charges of graft over black sand mining [5] suggest these operations persist due to corruption.Thus, ground monitoring of illegal activities is not only hard to apply at large scales but also dangerous for the involved individuals.Here we test whether remote sensing data can be used to identify and measure the scope of mining activities in the Philippines.
Black sand mining disturbs marine and coastal ecosystems and increases erosion and associated geohazards.Coastal black sand dunes serve as habitat for small animals and plants whose loss implies a threat to their predators.These dunes are also natural barriers to salt water, without them the sea often floods inland at high tides.The mining activities increase coastal erosion not only by directly removing sand but also by disrupting the sediment budget, often depriving areas in the down-current direction of their sand input.As a result, coastal erosion often continues to affect the areas even decades after cessation of the mining activities.This removal of material and associated erosion also likely results in land subsidence, which makes local communities particularly vulnerable to floods, damage from seasonal typhoons, and sea level rise.Additionally, black sand operations are often associated with groundwater extraction, which results in lowering of the water levels and increases the land subsidence.Subsidence is expected to continue after cessation of the water extraction due to residual compaction [6].Groundwater is used for processing (to wash sand, suppress dust, transport sand as a slurry, etc.) and for dewatering purposes for black sand to be extracted below the water table surface and for removing the water adhering to the sand grains [7].
We will use Interferometric Synthetic Aperture Radar (InSAR) to evaluate the vulnerability of local communities and the environment to black sand mining.InSAR has been used to successfully detect ground deformation associated with a number of geohazards, from earthquakes, to volcanic eruptions and land subsidence due to groundwater extraction [8][9][10][11][12][13][14][15].Here, we test whether the mining activities result in any detectable ground deformation despite their small spatial extent and evaluate for the first time the environmental effects of black sand mining with a high spatial resolution.To protect coastal ecosystems and communities, black sand mining is illegal in the Philippines within 200 meters of the shore [16] but the oversight of mining activities is poor.Using the high spatial resolution of the InSAR data we will evaluate the spatial impact of the mining activities.We will also characterize the magnitude of the associated subsidence, which will enable estimation of the time until the affected coastal areas are lowered to sea level elevation and exposed to flooding.
In this paper we aim to (1) demonstrate the value of remote sensing for monitoring of black sand mining activities and their environmental impact; (2) contribute a dataset of known and suspected mining areas in the Philippines; and (3) identify the areas of particular environmental and social concern subjected to increased hazards.

Identification of Mining Sites
The study of illegal activities presents a unique challenge since the behavior is intended to remain hidden.We use a multi-pronged approach to identify occurrences of black sand mining in the Philippines.First, we collected a list of known mining sites from various sources, including contacts on the ground, online news articles reporting either illegal mining activity or community protests related to illegal mining, and reports issued by the Philippine Mines and Geosciences Bureau (MGB) listing mining permits.Permits do not guarantee mining activity, but increase the likelihood by granting formal permissions to mining companies to work in the area.We also included potential sites based on iron ore or magnetite deposits reported in the United States Geological Survey (USGS) Mineral Resource Data System.
We then use optical satellite imagery to confirm the occurrence and precise locations of the mining sites, to characterize landscape changes over time, and to quantify the spatial extent of the mining activities.We rely on DigitalGlobe images provided by Google Earth and its historical imagery archive.These optical satellite data complement on-the-ground data by offering an objective form of measurement and enhanced spatial coverage of legal and illegal activities alike.

Evaluation of Environmental Effects
To identify the environmental impacts of black sand mining activities we rely on Interferometric Synthetic Aperture Radar (InSAR).InSAR allows us to characterize the extent and amplitude of land subsidence in areas with mining occurrence by enabling measurement of ground displacement in the radar line-of-sight (LOS) direction of a SAR satellite between different passes of the satellite over the same area [8].The Advanced Land Observing Satellite-1 (ALOS) of the Japanese Space Exploration Agency (JAXA) acquired SAR data with global coverage between late 2006 and mid 2011 on a 46-day repeat orbit, imaging most of the world's continents up to 25 times [17].We use ALOS data made available by the Alaska Satellite Facility (ASF) and process 14 ALOS frames on six tracks (Figure 1).Six SAR acquisitions are available on tracks 446, 447, and 448; nine acquisitions on track 449; and twenty-one acquisitions on tracks 445 and 450.
Remote Sens. 2016, 8, 100 3 of 16 imagery archive.These optical satellite data complement on-the-ground data by offering an objective form of measurement and enhanced spatial coverage of legal and illegal activities alike.

Evaluation of Environmental Effects
To identify the environmental impacts of black sand mining activities we rely on Interferometric Synthetic Aperture Radar (InSAR).InSAR allows us to characterize the extent and amplitude of land subsidence in areas with mining occurrence by enabling measurement of ground displacement in the radar line-of-sight (LOS) direction of a SAR satellite between different passes of the satellite over the same area [8].The Advanced Land Observing Satellite-1 (ALOS) of the Japanese Space Exploration Agency (JAXA) acquired SAR data with global coverage between late 2006 and mid 2011 on a 46-day repeat orbit, imaging most of the world's continents up to 25 times [17].We use ALOS data made available by the Alaska Satellite Facility (ASF) and process 14 ALOS frames on six tracks (Figure 1).Six SAR acquisitions are available on tracks 446, 447, and 448; nine acquisitions on track 449; and twenty-one acquisitions on tracks 445 and 450.  1) marked by white diamonds, and the location of the ALOS PALSAR data marked by the grey squares.The track numbers are specified below the respective tracks.Larges cities are indicated for reference with white circles.1) marked by white diamonds, and the location of the ALOS PALSAR data marked by the grey squares.The track numbers are specified below the respective tracks.Larges cities are indicated for reference with white circles.Table 1.Compilation of the sites of interest and their information.The first column shows the site identification, the second column its longitude and latitude, the third column whether the mining activities are legal (y) or illegal (n) or other type of activities, and the fourth column shows the source of information.The fifth column shows whether the InSAR mean velocity map has coherence in the area and the sixth column shows the detected LOS subsidence rates in cm/year.The last two columns show the name of the municipality with its population and the name of the province, respectively.The rows in grey correspond to sites out of the InSAR survey area or not associated with black sand mining activities.

Site
Lat First, ground displacement in the radar LOS direction is obtained from the phase difference of two SAR images acquired over the same area at different times (interferogram) [8].We use the ROI_PAC software [18] to produce 40 to 250 interferograms on each track.We remove topographic contributions using the Shuttle Radar Topography Mission (SRTM) 1-arc second digital elevation model [19].We co-register the interferograms of each frame to a master image, use the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU; [20]), and reference all interferograms to the same pixel.We subtract a ramp from each interferogram to remove any long wavelength signal (larger than 10 km, associated with orbital errors and interseismic deformation due to the subduction zone) so the obtained velocity maps highlight only localized deformation.
To precisely track ground deformation between the first and the last SAR acquisition we use the SBAS time series technique [21].The interferograms produced on each track are inverted to retrieve surface displacement through time relative to the first acquisition.We use only interferograms with small spatial and temporal baselines (i.e.position of the satellite close in space perpendicularly to the LOS at the different times (<1.6 km), and data acquired less than a year apart) to maintain high coherence (Figure 2).Since the PALSAR sensor onboard ALOS was L-band (wavelength of 23.6 cm), coherence is typically maintained longer than for X-and C-band sensors.Thus, we also add interferograms that span more than a year and cover the same season.The phase due to the digital elevation model (DEM) errors in SBAS time series is proportional to the perpendicular baseline history and, thus, it must be accounted for and corrected (Figure 2, right).We correct DEM errors in the in the time-domain following the method of Fattahi and Amelung, [22].Lastly, we perform the final pixel selection based on a temporal coherence threshold of 0.7 to eliminate pixels affected by phase-unwrapping errors [23][24][25].

of 17
First, ground displacement in the radar LOS direction is obtained from the phase difference of two SAR images acquired over the same area at different times (interferogram) [8].We use the ROI_PAC software [18] to produce 40 to 250 interferograms on each track.We remove topographic contributions using the Shuttle Radar Topography Mission (SRTM) 1-arc second digital elevation model [19].We co-register the interferograms of each frame to a master image, use the statistical-cost network-flow algorithm for phase unwrapping (SNAPHU; [20]), and reference all interferograms to the same pixel.We subtract a ramp from each interferogram to remove any long wavelength signal (larger than 10 km, associated with orbital errors and interseismic deformation due to the subduction zone) so the obtained velocity maps highlight only localized deformation.
To precisely track ground deformation between the first and the last SAR acquisition we use the SBAS time series technique [21].The interferograms produced on each track are inverted to retrieve surface displacement through time relative to the first acquisition.We use only interferograms with small spatial and temporal baselines (i.e.position of the satellite close in space perpendicularly to the LOS at the different times (<1.6 km), and data acquired less than a year apart) to maintain high coherence (Figure 2).Since the PALSAR sensor onboard ALOS was L-band (wavelength of 23.6 cm), coherence is typically maintained longer than for X-and C-band sensors.Thus, we also add interferograms that span more than a year and cover the same season.The phase due to the digital elevation model (DEM) errors in SBAS time series is proportional to the perpendicular baseline history and, thus, it must be accounted for and corrected (Figure 2, right).We correct DEM errors in the in the time-domain following the method of Fattahi and Amelung, [22].Lastly, we perform the final pixel selection based on a temporal coherence threshold of 0.7 to eliminate pixels affected by phase-unwrapping errors [23][24][25].The remaining signal contains noise contributions from atmospheric delay due to variation in water vapor content at different acquisition times.To evaluate the amplitude of the noise we follow the method of Chaussard et al. [12] and produce time series for pixels in non-deforming areas (metamorphic rocks).These time series provide a constraint on the amplitude of the atmospheric contamination at each epoch.Atmospheric noise can be as large as several centimeters in a single interferogram but is highly variable temporally and spatially, which implies that in non-deforming areas the trend through the time series has approximately a null slope.Thus, atmospheric noise corresponds to deviation from a longer trend, which represents deformation.Accordingly, we derive our estimates of ground subsidence using linear regressions through the time series, which are less affected by the atmospheric noise than single interferograms.
Chaussard et al. [13] showed that land subsidence due to groundwater pumping in Indonesia will result in flooding of urban areas within a few decades, this process being orders of magnitude The remaining signal contains noise contributions from atmospheric delay due to variation in water vapor content at different acquisition times.To evaluate the amplitude of the noise we follow the method of Chaussard et al. [12] and produce time series for pixels in non-deforming areas (metamorphic rocks).These time series provide a constraint on the amplitude of the atmospheric contamination at each epoch.Atmospheric noise can be as large as several centimeters in a single interferogram but is highly variable temporally and spatially, which implies that in non-deforming areas the trend through the time series has approximately a null slope.Thus, atmospheric noise corresponds to deviation from a longer trend, which represents deformation.Accordingly, we derive our estimates of ground subsidence using linear regressions through the time series, which are less affected by the atmospheric noise than single interferograms.
Chaussard et al. [13] showed that land subsidence due to groundwater pumping in Indonesia will result in flooding of urban areas within a few decades, this process being orders of magnitude larger than the rates of sea level rise.Here we calculate the approximate time for the mining sites areas associated with land subsidence to be at sea level elevation and, thus, permanently flooded and present this "time until flooding" so decision-makers can gain a sense of the time scales involved.This is the last step of our method, which is summarized in a flow chart in Figure 3.We evaluate the elevation based on the ASTER GDEM V2 Digital Elevation Model (DEM) and on the Global Multi-Resolution Topography (GMRT) [26].These 30 m DEMs lead to uncertainties of ~2 m on the relative elevation estimates for the mining sites.The time until flooding is equal to E/(S+R) with E the current ground elevation in cm, S the rate of subsidence in cm/year, and R is rate of sea level rise in cm/year.The two closest cities of Luzon Island that have current tide gauges monitoring (Legaspi and Davao) by the National Oceanic and Atmospheric Administration [27] show a mean sea level rise trend of 5.38 mm/year and 5.32 mm/year, respectively.These tide gauges are south of the region studied and can be affected by local subsidence e.g., [28].Thus, we verify that these sea level rise rates are in agreement with spatial altimetry data and non-Boussinesq models.Moon and Song [29] showed that between 1993-2010 Luzon Island was affected by a sea level rise of ~5 mm/year in agreement with the tide gauges data.Worldwide maps of sea level change based on satellite altimetry [30] show even faster sea level rise rates in the western Pacific at up to 10 mm/year east of Luzon island due to La Niña conditions (enhanced trade winds and cooling in the eastern tropical Pacific, corresponding to the absence of prolonged El Niño events in 1997-1998 and 2011).We, thus, consider that the tide gauge data are a good representation of the regional sea level rise for long time periods and use an average linear sea level rise of 5.35 mm/year in the Philippines.
Remote Sens. 2016, 8, 100 6 of 16 larger than the rates of sea level rise.Here we calculate the approximate time for the mining sites areas associated with land subsidence to be at sea level elevation and, thus, permanently flooded and present this "time until flooding" so decision-makers can gain a sense of the time scales involved.This is the last step of our method, which is summarized in a flow chart in Figure 3.We evaluate the elevation based on the ASTER GDEM V2 Digital Elevation Model (DEM) and on the Global Multi-Resolution Topography (GMRT) [26].These 30 m DEMs lead to uncertainties of ~2 m on the relative elevation estimates for the mining sites.The time until flooding is equal to E/(S+R) with E the current ground elevation in cm, S the rate of subsidence in cm/year, and R is rate of sea level rise in cm/year.The two closest cities of Luzon Island that have current tide gauges monitoring (Legaspi and Davao) by the National Oceanic and Atmospheric Administration [27] show a mean sea level rise trend of 5.38 mm/year and 5.32 mm/year, respectively.These tide gauges are south of the region studied and can be affected by local subsidence e.g., [28].Thus, we verify that these sea level rise rates are in agreement with spatial altimetry data and non-Boussinesq models.Moon and Song [29] showed that between 1993-2010 Luzon Island was affected by a sea level rise of ~5 mm/year in agreement with the tide gauges data.Worldwide maps of sea level change based on satellite altimetry [30] show even faster sea level rise rates in the western Pacific at up to 10 mm/year east of Luzon island due to La Niña conditions (enhanced trade winds and cooling in the eastern tropical Pacific, corresponding to the absence of prolonged El Niño events in 1997-1998 and 2011).We, thus, consider that the tide gauge data are a good representation of the regional sea level rise for long time periods and use an average linear sea level rise of 5.35 mm/year in the Philippines.We consider that the land subsidence rates observed are linear in time and persist even after the mining activities are over.This approximation is in part justified by the fact that coastal erosion often continues to affect the areas even decades after cessation of the activities due to sediment redistribution and by the fact that the land subsidence due to mining-related groundwater extraction is expected to continue due to residual compaction even after cessation of pumping [6].However, the linear extrapolation for time scales of half a century to a century may not be entirely appropriate as the mining sites are small and the aquifers involved are likely shallow and may consolidate leading to the subsidence to slow and stop.Thus, using the linear approximation may overestimate of the subsidence for long time scales and lead to lower end members of the values of time until flooding.We choose to give these lowest time values rounded to one standard deviation and present them as a "worst-case scenario" estimates representing what will happen if conditions are maintained as they are currently for decades.The development of new mining sites would result in new subsidence areas and could contribute in making these "worst-case scenario" estimates closer to being "average estimates" as sediment redistribution and aquifer compaction could become more wide-spread.

Identification of Mining Sites
We identified a total of twenty potential mining sites along the northern and northwestern coast of Luzon Island (Table 1).Precise coordinates for two mining sites were obtained from geo-tagged photos provided by contacts on the ground.Coordinates for six mining locations were identified from news articles referring to illegal black sand mining.Nine potential locations were identified from mining permits granted by the Mining and Geosciences Bureau in 2007, which were confirmed as mining sites with optical imagery.Locations of three coastal magnetite deposits were obtained from the USGS Mineral Resources Data System.However, for these magnetite deposits, the historical optical images did not reveal any clear mining activities for the 2003-2015 time span covered by the data available via Google Earth.Thus, out of the twenty potential sites, seventeen were confirmed as black sand mining sites with optical images (black labels in Table 1).
We use optical imagery to evaluate the amplitude of changes in landscape before and after mining is developed (Figure 4).We confirm that mining sites are very localized with perimeters on the order of a few hundreds of meters.Most of the sites are located on the coast, within the black sand dunes on the inland side of the beach (Figure 4).We also often notice the presence of water ponds nearby the extraction sites (Figure 4), which confirms the intensive use of water in the extraction process and the likely pumping of groundwater.
We consider that the land subsidence rates observed are linear in time and persist even after the mining activities are over.This approximation is in part justified by the fact that coastal erosion often continues to affect the areas even decades after cessation of the activities due to sediment redistribution and by the fact that the land subsidence due to mining-related groundwater extraction is expected to continue due to residual compaction even after cessation of pumping [6].However, the linear extrapolation for time scales of half a century to a century may not be entirely appropriate as the mining sites are small and the aquifers involved are likely shallow and may consolidate leading to the subsidence to slow and stop.Thus, using the linear approximation may overestimate of the subsidence for long time scales and lead to lower end members of the values of time until flooding.We choose to give these lowest time values rounded to one standard deviation and present them as a "worst-case scenario" estimates representing what will happen if conditions are maintained as they are currently for decades.The development of new mining sites would result in new subsidence areas and could contribute in making these "worst-case scenario" estimates closer to being "average estimates" as sediment redistribution and aquifer compaction could become more wide-spread.

Identification of Mining Sites
We identified a total of twenty potential mining sites along the northern and northwestern coast of Luzon Island (Table 1).Precise coordinates for two mining sites were obtained from geo-tagged photos provided by contacts on the ground.Coordinates for six mining locations were identified from news articles referring to illegal black sand mining.Nine potential locations were identified from mining permits granted by the Mining and Geosciences Bureau in 2007, which were confirmed as mining sites with optical imagery.Locations of three coastal magnetite deposits were obtained from the USGS Mineral Resources Data System.However, for these magnetite deposits, the historical optical images did not reveal any clear mining activities for the 2003-2015 time span covered by the data available via Google Earth.Thus, out of the twenty potential sites, seventeen were confirmed as black sand mining sites with optical images (black labels in Table 1).
We use optical imagery to evaluate the amplitude of changes in landscape before and after mining is developed (Figure 4).We confirm that mining sites are very localized with perimeters on the order of a few hundreds of meters.Most of the sites are located on the coast, within the black sand dunes on the inland side of the beach (Figure 4).We also often notice the presence of water ponds nearby the extraction sites (Figure 4), which confirms the intensive use of water in the extraction process and the likely pumping of groundwater.The main limiting factor for characterizing the development of new mining sites is the low temporal sampling of the optical images available through Google Earth.Nevertheless, our results The main limiting factor for characterizing the development of new mining sites is the low temporal sampling of the optical images available through Google Earth.Nevertheless, our results show that optical imagery is useful for characterizing legal and illegal mining activities, their scope, and their extent.With the increasing accessibility of high-resolution satellite data, monitoring of black sand mining activities from optical images could become a routine practice.In this study, we relied on manual evaluation of changes via visual inspection of optical images, and limited our study to cases with precise coordinates from our list of potential sites.In the future, the development of automatic classification using tools available through the Google Earth Engine could lead to automated and systematic evaluation of mining sites as soon as the optical imagery becomes available.

Evaluation of Environmental Effects
The mean InSAR velocity map produced (Figure 5) shows the speed at which the ground is moving with a high spatial resolution and for the entire coastal area of the Luzon Island.We produce a Google Earth product with the velocity maps to illustrate the ground deformation associated with black sand mining activities (Supplementary Material), which enables browsing of the large spatial area covered and easy sharing of our results with partners on the ground.show that optical imagery is useful for characterizing legal and illegal mining activities, their scope, and their extent.With the increasing accessibility of high-resolution satellite data, monitoring of black sand mining activities from optical images could become a routine practice.In this study, we relied on manual evaluation of changes via visual inspection of optical images, and limited our study to cases with precise coordinates from our list of potential sites.In the future, the development of automatic classification using tools available through the Google Earth Engine could lead to automated and systematic evaluation of mining sites as soon as the optical imagery becomes available.

Evaluation of Environmental Effects
The mean InSAR velocity map produced (Figure 5) shows the speed at which the ground is moving with a high spatial resolution and for the entire coastal area of the Luzon Island.We produce a Google Earth product with the velocity maps to illustrate the ground deformation associated with black sand mining activities (Supplementary Material), which enables browsing of the large spatial area covered and easy sharing of our results with partners on the ground.The background colors show the speed at which the ground is moving.Vertical motion corresponds to 1.2 times the LOS velocity.Blue colors show motion away from the satellite (subsidence) and red colors motion towards the satellite (uplift).White diamonds show the identified mining sites and red diamonds the subsiding sites detected in the velocity map but not originally identified (labeled InSAR sites).The insets show examples of subsiding areas (blue colors).White labels show sites with no coherence in the InSAR mean velocity map, cyan labels are sites with coherence but no detected subsidence, orange labels show subsiding mining sites with the subsidence rates in cm/year, and gray labels are areas with subsidence but not associated with black sand mining.
Decorrelation is the largest limitation of using InSAR for detection of subsidence associated with mining activities (Figure 5).Since the PALSAR sensor onboard ALOS was L-band (wavelength of 23.6 cm), it enables measurements through vegetation [12,15], but interferometry remains limited by temporal decorrelation of the radar signal; that is, changes in the characteristics of the ground between each satellite acquisition.When the ground characteristics change too much between satellites passes deformation cannot be evaluated.We found that seven out of the seventeen black sand mining sites surveyed with InSAR do not have sufficient coherence (Table 1, and Figure 5 sites marked as "no coherence").The majority of the mining sites are located in compressible alluviums, which outline most of Luzon Island and in coastal areas or near river beds with relatively flat topography and similar types of vegetation.Thus, differences in geology and terrain are unlikely to be the main contribution to the observation of decorrelation at some mining sites but not others.Additionally, some inland sites located in mountainous areas preserve correlation (e.g., site 4 and InSAR site 3) while coastal sites to the north lose correlation.The number of SAR acquisitions available likely influences the coherence, as the more frequent the acquisitions are, the more likely ground characteristics are to be maintained between each satellite pass.The ALOS satellite was acquiring on a 46-day orbit but acquisitions were not systematically collected worldwide.As a result, the number of acquisitions per track varies.The two tracks with the most SAR acquisitions are the tracks 445 and 450 (twenty-one dates) and track 449 (nine dates).Three tracks have a lower amount of SAR acquisitions: track 446, 447, and 448 (six dates).We notice that these three tracks with the lowest amount of SAR dates have a high number of sites suffering from decorrelation (three on track 448, one on track 447, and two on track 446) while tracks with higher numbers of SAR acquisitions (449 and 450) have no sites suffering from decorrelation.This suggests that with more frequent SAR acquisitions (with for example the 14-day repeat of the L-band ALOS-2 satellite, and the planned 12-day repeat of the NiSAR mission) more sites would have coherence and could be monitored for subsidence.The exception being track 445, which has a large amount of SAR acquisitions but where two mining sites suffer from decorrelation.It is worth noting that the overall correlation for track 445 is more than for tracks 446, 447, and 448, the majority of the area being correlated.The decorrelation at site 3 (small island north of Luzon) can be explained by the significant topography and vegetation, which differentiate it from most of the other coastal sites.One explanation for the occurrence of decorrelation in areas with large amount of SAR acquisitions (such as site 6 on track 445) would be either very large deformation or recent mining development but, unfortunately, we do not have access to detailed mining data on the timing of the activities and their amplitude to investigate which parameters control the loss of coherence.
For each site with coherence we calculate a mean subsidence rate, which corresponds to the slope of the best-fitting linear regression through the time series to limit the effect of atmospheric noise (Figure 6).The time series correspond to averages over the pixels experiencing subsidence.Out of the ten mining sites with coherence, four show no clear sign of deformation (from north to south sites 14, 11, 17, and 16) (Figures 5 and 6 first and second rows).Six out of the ten sites with coherence show subsidence at rates ranging from 1.3 cm/year to up to 4.6 cm/year in the Radar LOS (Table 1, Figures 5 and 6).Since InSAR is more sensitive to vertical than horizontal movements we convert LOS (dLOS) into vertical displacement (dv) for every time series using the average ALOS incidence angle (θ = 34.3˝): dv = dLOS/cosθ.Vertical ground displacement is 21% more than LOS displacement, i.e., 1 cm of LOS displacement corresponds to 1.2 cm of vertical displacement.The vertical subsidence associated with mining thus reaches up to 5.7 cm/year.
The InSAR velocity map also reveals subsidence in six additional sites (red diamonds on Figure 5) not originally identified.Out of these six sites, four were confirmed as mining sites with optical images (sites InSAR-1, 2, 3, and 5), but one was associated with mining activities other than black sand (gold extraction, site InSAR-2, located east of site 17) (Figure 5, Table 1).Subsidence in the radar LOS at the gold mining site reaches up to 4.5 cm/year (Figure 6; 5.4 cm/year vertically) and at the three additional black sand mining sites up to 2.6 cm/year (3.1 cm/year vertically) (Figure 6).These results show that InSAR is a good complement to other methods for identifying mining sites, but the distinction between black sand mining and other processes cannot be made based solely on the InSAR data (Figure 3).For example, based on optical images we confirmed that subsidence was also detected on volcanic deposits (InSAR site 4) and associated with agricultural activities (InSAR site 5 and Aringay).
Remote Sens. 2016, 8, 100 10 of 16 results show that InSAR is a good complement to other methods for identifying mining sites, but the distinction between black sand mining and other processes cannot be made based solely on the InSAR data (Figure 3).For example, based on optical images we confirmed that subsidence was also detected on volcanic deposits (InSAR site 4) and associated with agricultural activities (InSAR site 5 and Aringay).The other sites are ordered from north to south.The quadrants labeled "InSAR" are sites displaying subsidence but which were not originally identified as black sand mining sites.The last three InSAR sites are not associated with black sand mining activities (gray).
The InSAR mean velocity map provides both a high spatial coverage (entire northwest coast of Luzon Island) and a high spatial resolution, which allows us to evaluate the spatial extent of the subsiding areas compared to the mining sites extent characterized from optical images (Figure 7).We observe that the subsidence patterns extend much beyond the mining sites with an approximate ~10-to 100-times larger spatial extent (mining sites mostly extend on less that 0.05 km 2 while subsiding areas cover up to 5-10 km 2 , Figure 7).This larger extent of the subsiding areas likely reflects groundwater usage or sediment redistribution around the affected areas.We also notice that the spatial extent of the subsiding areas associated with the black sand mining sites is small compared to the extent of the land subsidence associated with other processes such as ground water extraction, e.g., [13][14][15]; oxidation of peatlands, e.g., [31]; natural sediment compaction e.g., [32]; slope motion e.g., [33]; sinkholes e.g., [34]; and hydrothermal or magmatic activity, e.g., [12,35].Thus, the small spatial extent of the land subsidence associated with black sand mining activities can be used to detect black sand mining sites and isolate this process from others.  5.Each dot marks a SAR acquisition time, the subsidence rates labeled correspond to the slope of the best fitting linear regression though each time series.The site names are labeled at the bottom left of each quadrant.The top four quadrants correspond to black sand mining sites with no clear deformation.The other sites are ordered from north to south.The quadrants labeled "InSAR" are sites displaying subsidence but which were not originally identified as black sand mining sites.The last three InSAR sites are not associated with black sand mining activities (gray).
The InSAR mean velocity map provides both a high spatial coverage (entire northwest coast of Luzon Island) and a high spatial resolution, which allows us to evaluate the spatial extent of the subsiding areas compared to the mining sites extent characterized from optical images (Figure 7).We observe that the subsidence patterns extend much beyond the mining sites with an approximate ~10-to 100-times larger spatial extent (mining sites mostly extend on less that 0.05 km 2 while subsiding areas cover up to 5-10 km 2 , Figure 7).This larger extent of the subsiding areas likely reflects groundwater usage or sediment redistribution around the affected areas.We also notice that the spatial extent of the subsiding areas associated with the black sand mining sites is small compared to the extent of the land subsidence associated with other processes such as ground water extraction, e.g., [13][14][15]; oxidation of peatlands, e.g., [31]; natural sediment compaction e.g., [32]; slope motion e.g., [33]; sinkholes e.g., [34]; and hydrothermal or magmatic activity, e.g., [12,35].Thus, the small spatial extent of the land subsidence associated with black sand mining activities can be used to detect black sand mining sites and isolate this process from others.We calculate, for the coastal sites with detected subsidence, the time until which the elevation will reach sea level and, thus, the area will be permanently flooded (Table 2).We follow the method described earlier and provide worse-case scenario estimates considering linear rates of land subsidence and of sea level rise (5.35 mm/year) to show what will happen if the conditions are maintained as they are currently for decades.We use the lowest bound of the current elevation provided by the DEM for these estimates and round them to one standard deviation.Sites with subsidence rates of 1.8 and 3 cm/year are projected to be underwater in ~50-70 years.Sites with subsidence rates of 4.3 and 4.6 cm/year are projected to be underwater in approximately 30-40 years.These results demonstrate the vulnerability of the coastal areas to black sand mining with only a few decades before coastal areas undergoing mining are permanently flooded.We calculate, for the coastal sites with detected subsidence, the time until which the elevation will reach sea level and, thus, the area will be permanently flooded (Table 2).We follow the method described earlier and provide worse-case scenario estimates considering linear rates of land subsidence and of sea level rise (5.35 mm/year) to show what will happen if the conditions are maintained as they are currently for decades.We use the lowest bound of the current elevation provided by the DEM for these estimates and round them to one standard deviation.Sites with subsidence rates of 1.8 and 3 cm/year are projected to be underwater in ~50-70 years.Sites with subsidence rates of 4.3 and 4.6 cm/year are projected to be underwater in approximately 30-40 years.These results demonstrate the vulnerability of the coastal areas to black sand mining with only a few decades before coastal areas undergoing mining are permanently flooded.

Discussion
Combining optical images and InSAR analysis allows for an in-depth identification of legal and illegal mining activities and their environmental impacts.To facilitate rapid response to illegal activities these data should be made accessible quickly after acquisition and processing should be automated.Implementation of systematic identification of new mining sites based on optical images is in the process of being developed in the Google Earth Engine.Currently, no SAR satellite provides global coverage of the Philippines.The Sentinel-1A, ALOS-2, and future NiSAR satellites should eventually allow for global Earth coverage once the missions are fully operational.The L-band sensor onboard the ALOS-2 and future NiSAR satellites would be particularly adapted for continuous monitoring of mining activities and their impacts.Thus, rapid, systematic processing of the SAR data that are planned to be distributed quickly after acquisition could lead to remote monitoring of illegal mining activities worldwide.
Our analysis reveals subsidence at several cm/year collocated with mining activities.Processes such as tectonic subsidence, natural compaction, and artificial settlement do not account for more than a few mm/year of subsidence (2-5 mm/year [32,36]).Additionally, natural subsidence occurs uniformly over large areas with consistent sediment deposits.In the Philippines InSAR reveals subsidence rates of ten times greater magnitude and with relatively small spatial extent.Thus, these subsidence rates cannot be explained by natural processes and likely result from man-made activities [13].Agricultural activities often contribute to rapid land subsidence [13,15], as confirmed at the InSAR site 5 and Aringay, labeled as fields and experiencing subsidence collocated with agricultural activities.However, with the exception of the subsiding the site 15, the mining sites are located close to the shore and not in the direct proximity of agricultural activities.Thus, groundwater extraction for agricultural activities are unlikely the cause of the observed subsidence.The colocation of the mining sites and the subsiding areas suggests a causal relationship with either the removal of material or associated effects, such as ground water pumping.
However, we observe subsidence associated with some mining sites but not all, which makes the causal relationship difficult to validate.Geology affects whether extraction results in land subsidence.For example, groundwater extraction results in compaction and subsidence only in compressible deposits [13,15].In the Philippines, we notice that the majority of the mining sites (subsiding and not) are located in compressible alluviums, which outline most of Luzon Island (beige in Figure 8).Thus, the difference between subsiding and non-subsiding sites cannot be explained by differences in geology.Two of the black sand mining sites of this study are not located on compressible deposits and they experience subsidence (red squares in Figure 8).It is likely that some surficial compressible deposits exist at these sites but have not been precisely mapped.
In mining areas where land subsidence was not detected it is possible that the rate of mineral and water extraction was not great enough to trigger subsidence above our detection threshold (cm/year).Recent data collection in the field confirms that a significant increase in illegal mining activity occurred between 2010 and 2013 (Table 3), unfortunately outside the time covered by the ALOS data.Continuous monitoring with new satellite platforms would help to better constrain the relationship between mining and land subsidence.To better constrain the parameters that influence the occurrence of land subsidence additional data are needed.In particular, the amount of subsidence likely depends on the time since the development of the mining activities, the amount of mining, and the amount of water extraction.Unfortunately, such data were not available to us due to the limited control on the mining activities.Future collaborations with communities on the ground could help access such datasets and further constrain the parameters controlling the amount and extent of land subsidence associated with mining.
The observation that mining is not always associated with subsidence does not allow us to validate the direct causality between mining and subsidence.We are however able to identify communities exposed to very high risk of flooding by typhoons and sea level rise.Escalation of legal or illegal mining activities in these areas may further exacerbate this risk.To better constrain the parameters that influence the occurrence of land subsidence additional data are needed.In particular, the amount of subsidence likely depends on the time since the development of the mining activities, the amount of mining, and the amount of water extraction.Unfortunately, such data were not available to us due to the limited control on the mining activities.Future collaborations with communities on the ground could help access such datasets and further constrain the parameters controlling the amount and extent of land subsidence associated with mining.
The observation that mining is not always associated with subsidence does not allow us to validate the direct causality between mining and subsidence.We are however able to identify communities exposed to very high risk of flooding by typhoons and sea level rise.Escalation of legal or illegal mining activities in these areas may further exacerbate this risk.

Conclusions
Remote sensing data are critical to monitor, control, and respond to black sand mining activities and their environmental and societal impacts.We show that optical images can be used to identify mining sites and respond to development of illegal activities.On Luzon Island we identified a total of twenty black sand mining sites, seventeen identified from contacts, news, and permits sources and confirmed with optical images and three identified in the InSAR mean velocity map and confirmed with optical images.We demonstrate that InSAR data can be used to characterize the land subsidence associated with the mining activities despite the small spatial extent of the mining sites.Seven of the twenty sites had no coherence in the InSAR data, thus no deformation could be measured, likely due to the limited number of acquisitions; four showed no ground subsidence above our detection threshold (cm/year); and nine experienced subsidence at rates ranging from 1.5 to 5.7 cm/year.The mean InSAR velocity map also highlights that the extent of the areas affected by subsidence is significantly larger than the mining sites themselves, likely associated with groundwater use or sediment redistribution around the affected areas.
Our results highlight the threat posed to coastal towns nearby black sand mining activities.Since most mining sites are at low elevation, the rapid subsidence results in high exposure to flooding and seasonal typhoons, and amplifies the effect of climate change-driven sea level rise.We show that several coastal areas will be at sea level elevation in a few decades due to the rapid subsidence.Since subsidence likely continues to affect the areas even decades after the cessation of mining activities due to the disruption of the sediment budget, characterization of the temporal evolution of land subsidence with longer SAR temporal coverage will be critical to mitigate environmental and societal effects of black sand mining activities.

Figure 1 .
Figure 1.Northern island of the Philippines (Luzon Island) with the identified sites of interest (Table1) marked by white diamonds, and the location of the ALOS PALSAR data marked by the grey squares.The track numbers are specified below the respective tracks.Larges cities are indicated for reference with white circles.

Figure 1 .
Figure 1.Northern island of the Philippines (Luzon Island) with the identified sites of interest (Table1) marked by white diamonds, and the location of the ALOS PALSAR data marked by the grey squares.The track numbers are specified below the respective tracks.Larges cities are indicated for reference with white circles.

Figure 2 .
Figure 2. Example for Track 445 of a perpendicular baselines-temporal baselines plot of a triangulated network of interferometric pairs.The dots show the SAR acquisitions dates.The left plot shows the interferograms produced on this track (black lines) after applying the baseline thresholds.The right plot shows the perpendicular baseline history of the set of SAR acquisitions, which is proportional to the DEM errors.

Figure 2 .
Figure 2. Example for Track 445 of a perpendicular baselines-temporal baselines plot of a triangulated network of interferometric pairs.The dots show the SAR acquisitions dates.The left plot shows the interferograms produced on this track (black lines) after applying the baseline thresholds.The right plot shows the perpendicular baseline history of the set of SAR acquisitions, which is proportional to the DEM errors.

Figure 3 .
Figure 3. Flow chart summarizing the methods used and results obtained in our survey to characterize black sand mining activities and their environmental impacts in the Philippines.Optical images and InSAR time series data are combined and show that remote sensing data provides an objective, reliable, safe, and cost-effective way to monitor black sand mining activities.Optical satellite data lead to identification of legal and illegal mining sites and InSAR allows us to evaluate the environmental impacts of black sand mining in terms of land subsidence.The InSAR time series method is described in detail in the method section.

Figure 3 .
Figure 3. Flow chart summarizing the methods used and results obtained in our survey to characterize black sand mining activities and their environmental impacts in the Philippines.Optical images and InSAR time series data are combined and show that remote sensing data provides an objective, reliable, safe, and cost-effective way to monitor black sand mining activities.Optical satellite data lead to identification of legal and illegal mining sites and InSAR allows us to evaluate the environmental impacts of black sand mining in terms of land subsidence.The InSAR time series method is described in detail in the method section.

Figure 4 .
Figure 4. Example of optical images illustrating the coastal landscape before (left) and after (right) the development of black sand mining activities (site 1).The mining activities are located within the black sand dunes, on the inland side of the beach.The white arrow highlights a water pond developed near the mining site.

Figure 4 .
Figure 4. Example of optical images illustrating the coastal landscape before (left) and after (right) the development of black sand mining activities (site 1).The mining activities are located within the black sand dunes, on the inland side of the beach.The white arrow highlights a water pond developed near the mining site.

Figure 5 .
Figure 5. InSAR mean LOS ground velocity map.The background colors show the speed at which the ground is moving.Vertical motion corresponds to 1.2 times the LOS velocity.Blue colors show motion away from the satellite (subsidence) and red colors motion towards the satellite (uplift).White diamonds show the identified mining sites and red diamonds the subsiding sites detected in the velocity map but not originally identified (labeled InSAR sites).The insets show examples of subsiding areas (blue colors).White labels show sites with no coherence in the InSAR mean velocity map, cyan labels are sites with coherence but no detected subsidence, orange labels show subsiding

Figure 5 .
Figure5.InSAR mean LOS ground velocity map.The background colors show the speed at which the ground is moving.Vertical motion corresponds to 1.2 times the LOS velocity.Blue colors show motion away from the satellite (subsidence) and red colors motion towards the satellite (uplift).White diamonds show the identified mining sites and red diamonds the subsiding sites detected in the velocity map but not originally identified (labeled InSAR sites).The insets show examples of subsiding areas (blue colors).White labels show sites with no coherence in the InSAR mean velocity map, cyan labels are sites with coherence but no detected subsidence, orange labels show subsiding mining sites with the subsidence rates in cm/year, and gray labels are areas with subsidence but not associated with black sand mining.

Figure 6 .
Figure 6.Ground displacement (in cm) as a function of time for the sites shown in Figure5.Each dot marks a SAR acquisition time, the subsidence rates labeled correspond to the slope of the best fitting linear regression though each time series.The site names are labeled at the bottom left of each quadrant.The top four quadrants correspond to black sand mining sites with no clear deformation.The other sites are ordered from north to south.The quadrants labeled "InSAR" are sites displaying subsidence but which were not originally identified as black sand mining sites.The last three InSAR sites are not associated with black sand mining activities (gray).

Figure 6 .
Figure 6.Ground displacement (in cm) as a function of time for the sites shown in Figure5.Each dot marks a SAR acquisition time, the subsidence rates labeled correspond to the slope of the best fitting linear regression though each time series.The site names are labeled at the bottom left of each quadrant.The top four quadrants correspond to black sand mining sites with no clear deformation.The other sites are ordered from north to south.The quadrants labeled "InSAR" are sites displaying subsidence but which were not originally identified as black sand mining sites.The last three InSAR sites are not associated with black sand mining activities (gray).

Figure 7 .
Figure 7. Example of an optical satellite image (a) showing the location and extent of the mining activities (black rectangles) and the corresponding InSAR mean velocity (b) for Site 4; (c) Mean times series for the maximum subsiding area and its best-fitting linear subsidence rate; (d) Example of profiles through the subsiding area showing the LOS velocity as a function of the distance along the profile (location shown with red dashed lines on (b).The top profile (1) is oriented NW-SE and the bottom profile (2) N-S.The black dots show the InSAR data, the gray area the standard deviation, and the red line the mean value.The horizontal black lines show the spatial extent of the mining site projected on the profile compared to the extent of the subsiding area.

Figure 7 .
Figure 7. Example of an optical satellite image (a) showing the location and extent of the mining activities (black rectangles) and the corresponding InSAR mean velocity (b) for Site 4; (c) Mean times series for the maximum subsiding area and its best-fitting linear subsidence rate; (d) Example of profiles through the subsiding area showing the LOS velocity as a function of the distance along the profile (location shown with red dashed lines on (b).The top profile (1) is oriented NW-SE and the bottom profile (2) N-S.The black dots show the InSAR data, the gray area the standard deviation, and the red line the mean value.The horizontal black lines show the spatial extent of the mining site projected on the profile compared to the extent of the subsiding area.

Figure 8 .
Figure 8. Map of the surficial geology of Luzon Island with the black sand mining sites identified with diamonds (geological data from the MGB of the Philippines).Labeled in white are mining sites with no coherence in the InSAR mean velocity.Labeled in cyan are the sites with coherence but no subsidence.Labeled in orange are the subsiding black sand mining sites with the rates indicated in cm/year.The two subsiding sites located in incompressible deposits are marked by red squares.

Figure 8 .
Figure 8. Map of the surficial geology of Luzon Island with the black sand mining sites identified with diamonds (geological data from the MGB of the Philippines).Labeled in white are mining sites with no coherence in the InSAR mean velocity.Labeled in cyan are the sites with coherence but no subsidence.Labeled in orange are the subsiding black sand mining sites with the rates indicated in cm/year.The two subsiding sites located in incompressible deposits are marked by red squares.

Table 2 .
Projected time in years until a site reaches sea level and is exposed to permanent flooding.Column 1 is the site name; Column 2 the LOS subsidence rate; Column 3 the subsidence rate converted to vertical; Column 4 is the vertical subsidence taking into account the mean predicted sea level rise; Column 5 is the mean elevation above sea level in meters; Column 6 is the estimated time until the site is at sea level elevation based only on subsidence rate; Column 7 is the time until the site is at sea level elevation considering subsidence and sea level rise.

Table 2 .
Projected time in years until a site reaches sea level and is exposed to permanent flooding.Column 1 is the site name; Column 2 the LOS subsidence rate; Column 3 the subsidence rate converted to vertical; Column 4 is the vertical subsidence taking into account the mean predicted sea level rise; Column 5 is the mean elevation above sea level in meters; Column 6 is the estimated time until the site is at sea level elevation based only on subsidence rate; Column 7 is the time until the site is at sea level elevation considering subsidence and sea level rise.

Table 3 .
Number of villages per province (from north to south) with reported illegal mining activity for different time periods.Due to timing and information constraints, we are unable to report mining activity for La Union and Pangasinan.