Land Surface Temperature Variation Following the 2017 Mw 7.3 Iran Earthquake

During an earthquake, crustal deformation, fluid flow, and temperature variation are coupled; however, earthquake-related land surface temperature (LST) variations remain unclear. To determine whether post-seismic fluid migration can cause changes in LST, and taking the Mw 7.3 2017 Iran earthquake as an example, we modeled surface cooling (CA) and warming (WA) areas induced by co-seismic slip and fluid migration using a thermo-hydro-mechanical (THM) coupled numerical simulation. Moreover, using nighttime LST data with 15-min resolution, the daily attenuation coefficient k of nighttime LST was extracted by attenuation function fitting, and the trend of the k time series was analyzed using the Mann–Kendall and Sen’s methods. Based on the comparison of k trends between the post-seismic and 2010–2016 periods, we obtained cooling and warming trends for the modeled CA and WA. The numerical simulation and observational data show good consistency, and both indicate that fluid migration caused by crustal deformation can lead to changes in LST. The numerical simulations show that after the Iran earthquake, the surface projection area of co-seismic slip correlated with a cooling area (CA), while the surrounding area correlated with a warming area (WA). For the LST observational data, the post-seismic k trends of the calculated CA and WA are positive and negative, indicating sustained cooling and warming processes, respectively. This study provides evidence that LST variation is caused by co-seismic crustal deformation and fluid migration and reveals the coupled evolution of deformation, fluid, and temperature fields. The results provide new insights into the mechanisms of seismic thermal anomalies.

To date, numerous observations of land surface temperature (LST) variation related to earthquakes have been reported [20][21][22][23], and fluid migration was considered as one of the possible contributors [20,24,25].However, debate remains as to the reliability of seismic related LST variations and their causes [20,22,25].Firstly, most seismic surface temperature research has focused on seismic thermal precursors with high uncertainties in timing, location, and intensity, all of which relate to poorly constrained pre-seismic processes, focal depths, magnitudes and structural settings.Secondly, satellite data with low time resolutions (e.g., four times per day) have been used; however, these data, such as Moderate Resolution Imaging Spectroradiometer (MODIS) LST data, are susceptible to environmental (e.g., local weather) and artificial disturbances [26,27].
Investigating LST variations caused by post-seismic fluid migration can minimize uncertainties in the timing and location of LST changes.Unlike the pre-seismic process, the co-seismic process causes observable fault slip and possible surface deformation [28,29]; this can lead to prolonged post-seismic fluid activity [9,10] and provides a long-term, high-intensity window for studying associated LST changes.Furthermore, areas of post-seismic LST variations caused by fluid migration can be modeled based on the thermo-hydro-mechanical (THM) coupled theory [30][31][32], which has been widely used in simulating nuclear waste storage, geothermal development, oil and gas exploitation, induced earthquakes, and earthquake fault slip behavior [17][18][19].
The use of high temporal resolution LST data can effectively improve the reliability of observed LST variations.The 15-min-resolution LST data, which contain the effects of solar forcing during the day and the decay process at night, are a valuable indicator of climate change [33,34], and contain more comprehensive temperature change information than data of lower time resolution [33,35].Extracting the daily characteristic thermal parameters of 15-min-resolution nighttime LST data has shown advantages in recognizing low frequency seismic-related signals and in minimizing the effects of high frequency disturbances [36].
To investigate whether fluid migration caused by crustal deformation can lead to changes in LST, we focused on the 2017 Mw 7.3 Iran earthquake, for which the 15-min-resolution LST data measured from Meteosat Second Generation (MSG) system are available [37].First, the thermo-hydro-mechanical (THM) coupled simulation method was used to calculate the distributions of cooling (CA) and warming (WA) areas caused by co-seismic slip and fluid migration.Then, we extracted the daily attenuation coefficient based on nighttime LST data, and compared the trends of daily attenuation coefficient time series for both CA and WA after the earthquake and during the period from 2010 to 2016.The study is one of the first attempts to investigate the influence of seismic related fluid migration on the LST by combining THM modeling and high temporal resolution LST analysis.Results of this study will improve our understanding of the coupling processes of earthquakes.

Setting and Data
On 12 November 2017, an Mw 7.3 earthquake occurred on the Iraq-Iran border [38] (hereafter referred to as the 2017 Iran earthquake; Figure 1).It was estimated that 600 to 700 people were killed, with more than 10,000 injured [39].After the earthquake, teleseismic and Interferometric Synthetic Aperture Radar (InSAR) data were used to investigate the co-seismic deformation and rupture processes, and an NE-dipping basement seismogenic fault was inferred [39,40].
Around the epicenter (white box area in Figure 1), prior to 2017 the most recent historical earthquakes (dated from 01 Jan 2009) included a series of Mw ≥ 5.0 events in 2013 among which the maximum is Mw 5.8; the 2017 earthquake was then followed by a series of Mw ≥ 5.0 aftershocks, including several Mw ≥ 6.0 events in 2018 and 2019 [38] (Figure 1; Table S1).Considering that the seismic moment of a Mw 6.0 event is ~0.01 times that of a Mw 7.3 event [41], LST data from 2014 to 2017 were used in this study; LST data from 2010 to 2013 were used for validation.
The LST data were downloaded from the Satellite Application Facility on Land Surface Analysis (LSA SAF, https://landsaf.ipma.pt/en/),where the retrieval of LST is based on clear-sky measurements from the MSG system in the thermal infrared window (MSG/SEVIRI channels IR10.8 and IR12.0) and with 3-km resolution at nadir.Theoretically, LST values can be determined 96 times per day from MSG [37].The Generalized Split-Window (GSW) algorithm [42] was chosen to retrieve LST.We used the MSG Toolbox Version 2.0 to process the data with sub-nadir pixel size and an error margin (uncertainty value) of less than 3.

THM Coupled Simulation Method
The THM coupled simulation method was used to calculate the thermal increment caused by fluid migration and the co-seismic slip of the fault so as to obtain warming and cooling areas of the surface.This coupling theory can be expressed by the hydro-mechanical (Equation ( 1)) and thermo-mechanical (Equation ( 2)) coupled equations [17,30,32,43]: In Equation ( 1), u is the Darcy velocity tensor, p f is the fluid pressure, ε vol is the volumetric strain of the porous medium, S is the storage coefficient, α B is the Biot-Willis coefficient, ∇ is the Hamiltonian operator, and ρ f is the fluid density and can be simplified if it is independent of spatial position.
In Equation ( 2), C p is the fluid heat capacity, q is the conductive heat flux, Q is the heat source or sink, and Q ted is the thermos-elastic damping.
Compared with the convective heat flux, the conduction term and the thermo-elastic damping term are far smaller, and Equation ( 2) can be simplified as: Here, this study focused on the distribution and relative strengths of the thermal increment; as such, the results shown in Sections 4.1 and 5.3 are normalized.Detail equations and parameters are shown in the Supplementary Materials (S1.1.and S1.2.; Table S2).

Attenuation Function Fitting of Nighttime LST Data
The attenuation function was used to fit of nighttime LST data to quantitatively evaluate the CA and WA obtained by numerical simulation.Based on the work of Göttsche and Olesen [44] and Inamdar et al. [33], we constructed two denudation functions for fitting the nighttime LST data: hyperbolic decay (f1; Equation ( 4)) and exponential decay (f2; Equation ( 5)): where t is time, T 0 is temperature around sunset, δT is the difference between T 0 and T(t) when t tends to positive infinity, t 0 is the time of sunset, and k is the attenuation coefficient in hours.The attenuation coefficient k can be used to indicate the strength of the cooling or warming action.Theoretically, an increase in k indicates an enhanced cooling effect, while a decrease indicates an enhanced warming effect; this was confirmed by applying f1 to observational data for 27 September 2017 (Figure 2).The results show that with enhanced cooling, k is obviously increased (Figure 2a), and with warming, k gradually decreases (Figure 2b).In the fitting process, t 0 is a known quantity, while T 0 , k, and δT are unknown variables.When there are data from within 2 h post-sunset, we considered the fitting to be effective, as our results show that missing data from within 2 h after sunset has little effect on the fit.In addition, to improve the fitting, we prescribed that the total number of data points used for fitting should not be less than 16.

Attenuation Coefficient Time Series Trends Analysis
Finally, trend analysis of the attenuation intensity time series within 1 month after the earthquake were conducted and compared the results with historical data to determine whether the study area has experienced cooling or warming.Specifically, we used the linear regression, Mann-Kendall trend test, and Sen's slope estimator to detect trends in the time series of attenuation coefficient k; these approaches have been widely used for the trend analysis of meteorological time series [45].

Mann-Kendall Trend Test
The Mann-Kendall test calculates the Z S of a time series; positive Z S indicates an increasing trend and vice versa [46,47], with Z S is given as: , S < 0 where where n is the number of data, sgn is the sign function, and where m is the number of tied groups and t i denotes the number of ties.A tied group is a set of sample data having the same value [45].

Sen's Slope Estimator
Sen [48] used the median of the set of slopes Q med to estimate the data trend, which is shown as: where the Q med sign reflects the trend and the value reflects the steepness of the trend.
In this study, we were interested in the differences between the trends of CA and WA, rather than the significance of the trend in any specific area.Therefore, we did not consider the significance level of the Mann-Kendall test or the confidence interval of Sen's slope estimator.Regardless, the methods and results of significance tests are presented in the Supplementary Materials (S2).

Co-Seismic Slip-Induced Surface Thermal Variation
Based on the slip model of Feng et al.
[49], we calculated the fluid convection induced thermal variation after the earthquake.The co-seismic slip of the 2017 Iran earthquake was mainly concentrated on two asperities with a maximum slip of 6 m (Figure 3a).In our THM model, co-seismic slip of the seismic fault can cause surface deformation, and elastic deformation can cause the convection of near-surface fluid, which can increase or decrease the heat at the surface.The modeled convection continued after the earthquake; specifically, cooling (i.e., CA) matches the surface projection of the slip area of the fault and two maximum cooling areas correspond to two asperities on the fault, while warming (i.e., WA) is observed around it (Figure 3b).We then selected one maximum cooling area (CA1) and a neighboring maximum warming areas (WA1) of the same size according to the THM model (Figure 1), to investigate the trace of convection-induced thermal variation using LST data.

LST Attention Intensity Trend Analysis
The f1 and f2 were used to fit the average LST in CA1 and WA1 from 12 November to 11 December for each year between 2014 and 2017 (Figure 4).It was found that the fitting of f1 was better than f2, which is consistent with the results of Inamdar et al. [33], and so trend analysis was performed for the k time series of f1.k trends detected by linear regression, Mann-Kendall trend test, and Sen's slope estimator were presented in Figure 5 and Table 1.k trends of both CA1 and WA1 are all negative except for that in CA1 following the earthquake in 2017.
For the linear regression, Mann-Kendall trend test, there were some differences in the k trends from 2014 to 2017.For example, from 2014 to 2015, the negative trend of the k value declined slightly, while that in 2016 was steeper.However, after the earthquake in 2017, the k trend in CA1 increased and was notably different from that in the historical period; in contrast, the trend in WA1 was within historical fluctuations.
If they were not affected by local cooling and warming, the adjacent CA1 and WA1 should have exhibited consistent k trends, much like the near-uniform trends seen for the period 2014-2016.However, after the 2017 Iran earthquake there was a notable difference between the trends for CA1 and WA1.The k trend in the CA1 region increased slightly, indicating an increase in the local cooling effect, while the k value in the WA1 region decreased steeply, indicating an opposite local thermal effect, i.e., the warming effect (Figure 5; Table 1); this shows good consistency with the numerical simulation results.

Reliability Analysis of k Trends
The comparison of nighttime LST attenuation trends in CA1 and WA1 for the post-seismic period shows that the CA1 and WA1 k trends only deviated following the 2017 Iran earthquake; this indicates cooling and warming processes, respectively, and is consistent with the numerical simulation results.Especially for CA1, the post-seismic attenuation coefficient k trend was remarkably larger than that for 2014-2016.We observed similar results for CA2, another cooling area adjacent to CA1 (Figure 1; Table S3).
The nighttime LST attenuation trends in CA1 and WA1 during the period from 2014 to 2016 shows that the CA1 and WA1 exhibited consistent k trends; this is consistent with the assumption that, for two neighboring areas, if there were no local cooling or warming processes, their k trends would be similar.In addition, the k time series during 2010-2013 shows similar characteristics for CA1 and WA1 (Table S3; Figures S1 and S2).And two other warming areas, i.e., WA2 and WA3, also show above characteristics (Figure 1; Table S3).
Furthermore, for the three warming areas, it is found that regression coefficients in 2017 might be related to the intensities of the warming increment.Specifically, the warming increment in WA1 is greater than that in WA3, in turn, greater than that in WA2.And the same is true for the regression coefficients of the three.Therefore, the deviation between regression coefficients of CA1 and WA1 following the 2017 earthquake indicate the increasing difference in LST of the two areas (Figure 5; Table 1).

Seismic Thermal Anomaly Mechanism
Owing to insufficient observational data and a poor understanding of pre-seismic processes, to date there has been no unified understanding of the mechanisms of pre-seismic thermal anomalies [26,27].Pre-seismic thermal anomalies were reported after the 2017 Iran earthquake [50,51].Akhoondzadeh et al. [50] identified a warming area of about ~3-degree size around the epicenter preceding the earthquake by 15 days, but both the pre-seismic anomalous area and duration were different from the post-seismic LST variation.However, the results of this study show that co-seismic crustal deformation can cause post-seismic fluid migration and LST changes, similar to the 'earth degassing' mechanism [24,25], a process that may be important for pre-seismic thermal anomalies.

Effect of Afterslip on LST Change
In addition to poro-elastic rebound, afterslip is an important mechanism affecting post-seismic deformation [52].Following the 2017 Iran earthquake, Feng et al. [49] revealed the surface deformation caused by afterslip of the seismic fault using InSAR.Their results showed that afterslip was mainly distributed on the up-dip fault plane adjacent to the co-seismic slip.The maximum cumulative slip 1 month after the earthquake was ~0.6 m (Figure S3a).We calculated the surface thermal increment caused by cumulative post-seismic slip and found that, under the same geological conditions, the thermal variation caused by afterslip has a wider range, but that the intensity did not exceed 0.1 times the co-seismic slip (Figure S3b).Therefore, although the co-and post-seismic temperature change areas are partially superimposed, we believe that the main contribution was from the co-seismic slip.

Limitations
It was a one-month period in which we conducted the trend analysis in this study.This arrangement is mainly made with the following concerns: sufficient data points are necessary for trend analysis as insufficient data could lead to abnormal results (Tables S3 and S4); post-earthquake deformation in the one month is adequately constrained [49]; and the warming and cooling effects are still considerable (Figure S4).It is still a preliminary study on the relation between the post-seismic fluid migration and LST changes.And improvement of this method or a more sufficient and robust method to detect the LST changes still need further works.
LST data in November 12 and 13, 2017, when the warming or cooling effect induced by the fluid migration have the maximum intensity based on the THM model (Figure S4), were not used in this study.Because their RMSEs and residual errors of the fitting are far greater than the other data, of which most RMSE are less than 0.5 and residual errors less than 1 K.We presented the comparison of fitting results in Figure S5 and S6.Besides, the LST data in the two days became discrete shortly after the main shock, we cannot distinguish if it is affected by other warming or cooling factors.However, based on the THM model (Figure S4), it can be inferred that lack of the first two days following the 2017 Iran earthquake have limited influence on the one-month-long trend analysis.

Conclusions
To determine whether post-seismic fluid migration can cause changes in LST, and taking the Mw 7.3 2017 Iran earthquake as an example, we modeled surface cooling (CA) and warming (WA) areas induced by co-seismic slip and fluid migration.Moreover, using nighttime LST data with 15-min resolution, the daily attenuation intensity coefficient k of nighttime LST was extracted by attenuation function fitting, and the trend of the k time series was analyzed using the Mann-Kendall and Sen's methods.Different k trends between CA and WA only occurred after the earthquake; that is, our results provide evidence that LST variation is caused by co-seismic crustal deformation and fluid migration.
THM numerical simulation reveals areas of surface cooling and warming near the epicenter of the 2017 Mw 7.3 Iran earthquake; these areas can be attributed to co-seismic slip and fluid migration.Following the 2017 Iran earthquake, the cooling area was mainly located along the surface projection of co-seismic slip along the fault, while the warming area surrounded the cooling area.Simulation results suggest that afterslip has little effect on this process; the one-month cumulative slip contributed less than 10% to the surface temperature variation.
The post-seismic cooling and warming areas were confirmed by analysis of nighttime LST data.The consistency between the numerical simulation and the observed data indicate that fluid migration induced by crustal deformation can cause changes in surface temperature.
The results of this study show that attenuation function fitting can be used to successfully extract temperature variation trends from data with a high temporal resolution.Such data contain obvious responses to warming and cooling processes, and can resist high frequency disturbances.A hyperbolic function can better fit the attenuation process of nighttime temperature as compared with an exponential function.Following the 2017 Iran earthquake, the attenuation intensity of the cooling area increased remarkably compared with the historical level.

Figure 1 .
Figure 1.(a) Regional location of the study area (inset map).(b) Tectonic setting of the study area, with topography from Shuttle Radar Topography Mission (SRTM) 90-m digital elevation model (DEM) data (http://srtm.csi.cgiar.org).(c) Earthquake time sequence recorded by the United States Geological Survey (USGS [38]) from 2009 to September 2019 in the white box.Red stars denote the epicenters, and white-red balls represent proposed focal mechanisms for the 2017 Iran earthquake from the USGS and the Global Centroid-Moment-Tensor (CMT).Red lines represent faults, with the thick red line representing the Zagros thrust belt.The central color map shows post-seismic temperature increments following the 2017 Iran earthquake obtained by numerical modeling; the white box shows the area of earthquake events; red colors indicate warming and blue colors indicate cooling.Dashed boxes CA1 and CA2 denote two cooling areas; WA1, WA2 and WA3 denote three warming areas.

Figure 2 .
Figure 2. Response of attenuation coefficient k to (a) cooling and (b) warming.The horizontal axis is time from sunset to sunrise the next day; the vertical axis is temperature.Black points and lines represent observed raw data and the best fit of the hyperbolic decay function, respectively.In (a), blue points and lines represent data points and fitting after superimposing the linear cooling of different intensities; red lines in (b) show the same for warming.A value of '−0.04 K/h' represent a temperature decreases 0.04 K per hour, and R is the root mean square error (RMSE).

Figure 3 .
Figure 3. (a) Co-seismic slip on the fault of the 2017 Iran earthquake and the geometric model.The red star denotes the USGS source position on the fault [38].(b) Post-seismic thermal increment at the surface induced by co-seismic slip.Cool colors indicate cooling and warm colors indicate warming; the red star denotes the projection of USGS hypocenter on the surface.

Figure 4 .
Figure 4. Nighttime land surface temperature (LST) data fitting following the 2017 Iran earthquake.Red plus signs show original data; (a) blue lines are the hyperbolic decay f1 fitting for CA1, and (b) red lines are the f1 fitting result for WA1.Black lines denote exponential decay f2 fitting, and R is the root mean square error (RMSE).

Figure 5 .
Figure 5. Historical k time series trends during (a) 2014, (b) 2015, (c) 2016; and (d) the post-seismic k time series trend in 2017.Blue triangles denote k in the cooling area (CA1) and blue lines denote the linear regression through the data.Red triangles denote k in the warming area (WA1) and red lines denote the linear regression through the data.Regression coefficients are given in unit of hours per day (h/d).MK, Mann-Kendall.
Figure S1: Historical k time series trends during (a) 2010, (b) 2011, (c) 2012; (d) 2013.Blue triangles denote k in the cooling area (CA1) and blue lines denote the linear regression through the data.Red triangles denote k in the warming area (WA1) and red lines denote the linear regression through the data.Regression coefficients are given in unit of hours per day (h/d).MK, Mann-Kendall.
Figure S2: Regression coefficients and confidence interval of k time series during 2010 to 2017.Red diamonds and error bars indicate the regression coefficients and confidence intervals of WA1 warming zone respectively, and blue ones indicate those of CA1. Figure S3: (a) Afterslip on the fault of the 2017 Iran earthquake and the geometric model.Red star denotes the United States Geological Survey (USGS) source position on the fault [38].(b) Thermal increment at the surface caused by afterslip.Warm colors indicate warming and cool colors indicate cooling; the red star denotes the USGS epicenter on the surface.
Figure S4: Evolution of post-seismic thermal increment at the center of WA1.Red circles represent the time steps of calculation.0 represents the origin time of the 2017 Iran earthquake. Figure S5: comparison of fitting results of CA1 in 2017.
Figure S6: comparison of fitting results of WA1 in 2017.Author Contributions: Data curation, Z.J.; formal analysis, C.Z., G.Z. and Y.L.; methodology, C.Z. and Z.J.; supervision, X.S.; visualization, Y.L.; writing -original draft, C.Z.; writing -review and editing, C.Z., X.S., G.Z. and Y.L. Funding: This research was funded by the National Natural Science Foundation of China (grant number 41631073), the National Key Research and Development Program of China (grant number 2018YFC1503602), and the China Seismic Experimental Site Project (grant number 2018CSES0205).

Table 1 .
Historical and post-seismic time series trend 1 .