Constraining Historical Earthquake Sequences With Coulomb Stress Models: An Example From Western Türkiye

Knowledge of earthquake source faults is crucial for the calculation of robust Coulomb stress models. However, source faults are often poorly constrained, especially for pre‐instrumental events, and these historical earthquakes are commonly studied with little or no consideration of other nearby events. We introduce an approach using Coulomb Stress Transfer (CST) modeling to investigate historical earthquake sequences and constrain possible source faults. Using historical and instrumental records from the Büyük Menderes Graben, western Türkiye, we create an ensemble of earthquake sequences featuring multiple rupture scenarios for individual earthquakes, and model both coseismic and interseismic CST. We filter and evaluate the models based on criteria to gain knowledge on historical earthquakes and their source faults and assess the current stress state and related seismic hazard of the investigated fault network. For our study area, the results provide further constraints on the source faults of several historical earthquakes, including the destructive MW 7.0 earthquake in 1899. The approach presented herein is applicable to other tectonically active areas where the causative faults of historical earthquakes are poorly constrained from existing data sets (e.g., paleoseismology, damage records), providing a new tool to help decipher historical earthquake sequences and improve modeling studies.

remains a major challenge.Other studies have addressed the issue of source fault uncertainty by modeling multiple possible source faults for individual earthquakes (e.g., León-Loya et al., 2022;Mohammadi et al., 2019) and assessing the differences in the resultant stress patterns.However, this approach does not assess the likelihood of each source fault scenario and would only work for a small number of earthquakes with multiple possible scenarios.
We develop an approach using Coulomb stress modeling to help constrain the source fault(s) for a sequence of historical earthquakes.We utilize the hypothesis that earthquakes promote (or delay) ruptures on nearby faults through CST (e.g., Harris & Simpson, 1992;Scholz, 2010;Stein, 1999;Stein et al., 1997).Instead of focusing on solely single earthquake(s), we model CST throughout an earthquake sequence and evaluate different patterns of stress transfers for a range of possible source faults.We propose that when a historical rupture is poorly constrained, we can use CST models of prior and subsequent earthquakes to constrain which source fault/segment is most likely to have ruptured.
We study an earthquake sequence in the Büyük Menderes Graben, located in the western Anatolian extensional province in SW Türkiye (Figure 1a).Western Türkiye and parts of the Aegean Sea are undergoing N-S extension, mostly accommodated across a series of ∼E-W trending normal faults/grabens, at a rate of ∼20 mm/yr (Aktug et al., 2009;Kurt et al., 2023;McClusky et al., 2000), driven by the westward escape of the Anatolian Microplate and roll-back at the Hellenic Arc subduction zone (Barka & Reilinger, 1997;Bozkurt, 2001;Jackson, 1994;McKenzie, 1972McKenzie, , 1978;;Taymaz et al., 2007).The Büyük Menderes Graben is >120 km long and has a record of at least five large earthquakes (M W > 5.8) over the past 400 years.It is therefore one of the largest and most active grabens in the region.The Büyük Menderes Graben fault (BMGF), an E-W trending, S-dipping normal fault zone forms the northern margin of the graben (Figure 1b).One of the largest and most destructive events in the region likely occurred on the BMGF in 1899, with an estimated M W 7.0 (Altunel, 1999; N. N. Ambraseys & Jackson, 1998;Stucchi et al., 2013), causing widespread destruction in the Büyük Menderes Graben and the neighboring Denizli Basin.Despite the availability of detailed historical sources, there is debate over which section of the BMGF ruptured in this earthquake.Several studies propose various locations, lengths and throws of the co-seismic fault scarp (Altunel, 1999; N. N. Ambraseys & Jackson, 1998;Ergin et al., 1967;Ocakoğlu et al., 2013) (see Supporting Information S1 for a summary).Given the extent of damaged areas and the magnitude of the earthquake, it is likely that the event ruptured at least one or two segments of the BMGF (Figure 1c).An alternative scenario proposed by Ocakoğlu et al. (2013) features two main shocks, one in the Büyük Menderes Graben and a second in the neighboring Denizli Basin.If they occurred within a short time frame of each other, it is unlikely that the two events would be distinguishable in the historical records.Since the historical record leaves several unanswered questions, we focus our study on this earthquake by studying the historical earthquake sequence using Coulomb stress modeling.
This study aims to; (a) present an approach to assess historical earthquakes and their source faults using Coulomb stress modeling, (b) provide new constraints on historical earthquakes in the Büyük Menderes Graben region, especially the debated 1899 earthquake, and (c) create a Coulomb stress model of the region to assess the current stress state of the fault network.

Active Faults and Earthquakes
The studied fault network is based on existing fault catalogs (R. Emre & Sözbilir, 2007;Ö. Emre et al., 2018;Koçyiğit, 2005), with some additional fieldwork observations.All the studied faults consist of multiple faults, segments, or branches, we have simplified them to single traces for modeling (Figure 1b and Table 1).
The evolution of the Küçük Menderes Graben is linked to uplift of the Bozdağ Horst between the S-dipping BMGF in the south and the N-dipping Gediz Graben fault in the north (Seyitoğlu & Işik, 2009).Like the other grabens, normal faults on both sides separate the graben fill from the metamorphic rocks of the Menderes Massif; however, active faults in the graben are poorly mapped.Tectonic activity is not proven, but large earthquakes with M W 6.2 and 6.3 have been recorded in the area in 1880 (Stucchi et al., 2013) and 1928 (Papazachos & Papazachou, 1997).The geomorphology of the Küçük Menderes Graben suggests an active fault zone along its southern margin, the easternmost part of which was mapped by R. Emre and Sözbilir (2007).We discovered small but well-preserved fault planes in the central part of the morphologic scarp, indicating normal faulting with a sinistral component (Figure 2).
The Denizli Basin is a NW-SE-trending complex graben structure hosting several active fault zones.The north-eastern side of the basin is formed by the Pamukkale fault, which caused destructive historical earthquakes, for example, in 60 AD (Hancock & Altunel, 1997;Hancock et al., 2000;Kumsar et al., 2016;Piccardi, 2007) and controls the activity of thermal springs (Brogi et al., 2014(Brogi et al., , 2016)).The south-western side of the graben is the Babadağ fault, while the graben center hosts the Kaleköy fault zone, consisting of multiple faults (names and exact locations vary among publications) and bounding the Cürüksu sub-graben in the center of the Denizli Basin (Kaymakci, 2006;Koçyiğit, 2005).We assume that the basin-bounding faults, the Pamukkale and Babadağ faults, are the dominant structures and the Kaleköy fault is limited by these.The last active fault incorporated into the model is the Honaz fault, which delimits the southern side of the Denizli Basin and ruptured in at least one ∼M W 5.5 earthquake in 1965 (Koçyiğit, 2005;Kumsar et al., 2016).Based on historical and recorded seismicity, geomorphology, and paleostress analyses, all these structures are active normal faults, capable of rupturing moderate to large earthquakes (Brogi et al., 2014;Irmak, 2013;Kaymakci, 2006;Koçyiğit, 2005;Özkaymak, 2015).
Although Anatolia has been settled for millennia and earthquake records date back thousands of years, the record is incomplete for most of this time.We study 16 earthquakes with M W > 5.5 (Supporting Information S1 and Table 2) that have occurred since the mid-17th century.Given the good historical records and dense population of the region, we assume that the earthquake record for this period, study area, and magnitude is complete.
We modeled the following earthquakes with one to three scenarios (S1-S3).A detailed description of each event, including references and explanations for the choice of scenarios, is given in Supporting Information S1.
A M W 6.6 earthquake occurred in the central Büyük Menderes Graben (BMG) in 1645 or 1646, which most likely ruptured the Atça segment of the BMGF (S1).A possible alternative would be the Umurlu segment (S2) or the Pamukören segment (S3) (N.Ambraseys, 2009;Kumsar et al., 2016;Papazachos & Papazachou, 1997).It was followed by a M W 6.7 shock in 1651 in the Denizli basin.According to Koçyiğit (2005) it occurred on the Pamukkale fault (S1), but evidence for this is sparse (N.Ambraseys, 2009;Papazachos & Papazachou, 1997), hence we model alternative scenarios on the Kaleköy fault (S2) and Babadağ fault (S3).In 1653 another earthquake of M W 6.7 ruptured the western part of the BMG (N.Ambraseys, 2009;Ocakoğlu et al., 2013;Papazachos & Papazachou, 1997;Stucchi et al., 2013).We model three scenarios, S1 rupturing solely the İncirliova segment and S2 and S3 rupturing different proportions of the İncirliova and the adjacent Umurlu segment.Two well documented large earthquakes occurred in the Denizli basin in 1702 and 1717 (N.Ambraseys, 2009;Koçyiğit, 2005;Kumsar et al., 2016;Ocakoğlu et al., 2013;Papazachos & Papazachou, 1997;Stucchi et al., 2013).The first one most likely ruptured the Kaleköy fault zone (S1) along its entire length, the second one likely the Babadağ (S1) or possibly the Pamukkale fault (S2).The next large earthquake is reported from the Küçük Menderes Graben in 1880 (Stucchi et al., 2013).We model it on the Küçük Menderes Graben Fault (KMGF) (S1), but due to the poor constraints on this event we include a second scenario where it does not occur on the modeled fault network at all (i.e., leaving it out).Two large earthquakes occurred near the Aegean coast in 1890 and 1891, likely at the Kuşadası fault zone (N.Ambraseys, 2009;Stucchi et al., 2013).We model one scenario each, the first event on the western segment and the second event in the eastern part.This was closely followed by a M W 6.6 earthquake in 1895 in the adjacent, western BMG, most likely on the Umurlu segment (S1), or possibly on the İncirliova segment (S2) or, less likely, parts of both segments (S3) (N.Ambraseys, 2009;Ocakoğlu et al., 2013;Papazachos & Papazachou, 1997;Stucchi et al., 2013).In 1897 a M W 6.2 earthquake occurred at the transition between the eastern BMG and the Denizli basin (Stucchi et al., 2013).Possible source faults are the eastern BMG (Pamukören segment, S1) or the Kaleköy fault (S2).This was followed by the devastating M W 6.9-7.0 1899 earthquake.As  previously discussed, data on the extent and location of surface ruptures vary a lot, but the earthquake likely ruptured more than one segment of the BMGF (N.Ambraseys, 2009;Ergin et al., 1967;Ocakoğlu et al., 2013;Papazachos & Papazachou, 1997;Stucchi et al., 2013).Our model includes two scenarios, rupturing the Atça segment and parts of the Pamukören segment (S1) or the entire Pamukören and Buharkent segments (S2).A M W 5.8 earthquake occurred in the same region in 1920 (KOERI, 2020), possibly rupturing the Atça (S1) or the    A final earthquake occurred in the BMG in 1986 (Leptokaropoulos et al., 2013;Taymaz, 1993), however due to the comparably small magnitude and poor constraints on its location we decided to exclude this scenario.

Defining Faults, Earthquakes, and Plausible Rupture Scenarios
We compiled data on the location, geometry, kinematics, and activity of active faults from literature and field observations (Figures 1b and 2, Table 1).The data consists of instrumental and pre-instrumental earthquakes of M W > 5.5 from earthquake catalogs and historical studies (Supporting Information S1), containing magnitudes, epicenter locations, macroseismic shaking information and other details.For each historical earthquake, one to three plausible rupture scenarios are chosen based on historical shaking, fault segmentation and expected rupture length (Figure 1c and Table 2 and Supporting Information S1).Scenarios may be different source faults/ segments, rupture parameters or non-occurrence.Earthquakes with reliable historical data are represented by a single scenario, whereas more uncertain events are represented by two or three scenarios.We model a maximum of three scenarios per earthquake to limit the number of possible combinations.

Coulomb Stress Modeling
We built the 3D fault network using the "3D-Faults" software (version 2.5), an improved version of the code published by Mildon et al. (2016).This version features more fault input formats, a user interface, and the option to detect and cut faults intersecting at depth.Faults with multiple segments (e.g., BMGF) are simplified into a single continuous fault, earthquakes can then be modeled on individual segments by specifying the position of the rupture segment along the fault zone (Figure 3b, Table 2).
Table 1 shows the geometry and slip rates used for modeling the fault zones with the 3D-Faults software.Where multiple values are available in the literature, we take the average.Faults are extended straight to a seismogenic depth of 18 km using the given dip angles, except for the BMGF, where a listric geometry was used, with dip angle varying between 65° (surface) and 23° (bottom), consistent with the general listric geometry of the large, regional detachment faults (Sengör et al., 1985;Çiftçi & Bozkurt, 2009;Çifçi et al., 2010).Faults were built with a grid size of 1 km, and where different faults intersect at depth, we assume that the largest graben bounding faults are the dominant structures and therefore we cut the other fault at the depth of intersection.The modeled faults consist of >8,000 rectangular elements, creating strike-variable fault surfaces with either listric or planar geometry to depth.For the Küçük Menderes Graben fault there are no published slip rates, therefore we assign a rate of 0.7 mm/yr, which is in the range of other faults in the region with similar dimensions and morphologic expression.
A concentric slip distribution is assigned to the source fault/segment with the maximum slip in the center of the rupture.Similar to other studies (e.g., Mildon et al., 2022), slip at the fault edges and base is set to zero, while slip at the surface is 40% of the maximum slip (Manighetti et al., 2007).Rupture length is inferred from magnitude-length relationships (Wells & Coppersmith, 1994) and the mapped length of the source fault/segment.We iterate the value of maximum slip at the center of the concentric slip distributions, so the calculated magnitude matches the catalog magnitude (Table 2).For large earthquakes, the entire seismogenic depth ruptures, otherwise if the rupture length is shorter than the down-dip-length the rupture area is assumed to have an aspect ratio of 1.We model a seismogenic thickness of 18 km, which is the maximum depth of most recorded seismicity in the region (Leptokaropoulos et al., 2013).We use a coefficient of friction of 0.4 which is appropriate for modeling tectonic earthquakes (e.g., King et al., 1994).'3D-Faults' output is compatible with 'Coulomb 3.4' (Toda et al., 2011), which we use to calculate coseismic CST for each rupture scenario.
We model coseismic CST for each of the specified rupture scenarios, resolving the CST onto the mapped fault surfaces and for receiver faults at a specified depth and geometry.We specify a receiver fault geometry similar to the BMGF of 180°/60°/−90° (dip direction/dip/rake), because most earthquakes/aftershocks triggered by CST occur on planes similarly oriented as the source fault (Steacy et al., 2005).We calculate the CST at 8 km depth, the approximate depth of most hypocenters in the region (Leptokaropoulos et al., 2013).
Furthermore, we include annual interseismic loading since 1600 which is calculated using the 'back-slip' method (Deng & Sykes, 1997;Savage, 1983;Verdecchia & Carena, 2015).This approach models 'virtual negative displacements', that is, incremental (millimeter sized) reverse slip applied on the normal faults, and then the CST is resolved for normal motion on the same fault surfaces.The slip increments are based on fault slip rates from the literature (Table 1).We assume that the slip rate(s) occur in the center of the fault (as few slip rates are given for specified locations) and interpolate a triangular slip rate profile along the fault, decreasing to zero at the ends.
Finally, we calculate the cumulative CST over time for each of the 5,184 possible earthquake sequences created by combining all rupture scenarios (Figure 3a).For the CST at a specified depth, we sum the coseismic CST only (Figure 3c).For the CST resolved onto mapped fault surfaces, we sum the coseismic and interseismic CST, assuming that stress on a ruptured segment decreases to zero after a rupture (e.g., Mildon et al., 2019;Verdecchia andCarena, 2015, 2016).These cumulative CST models are then used to perform the filtering and evaluation steps to constrain which rupture scenarios are most plausible.An example for a single modelled earthquake sequence is visualised in Movie S1.

Filtering and Evaluation
We aim to determine the models that are more plausible by applying assumptions about CST and fault interaction to filter out less-plausible rupture scenarios.
Studies of earthquake sequences have shown that subsequent earthquakes predominantly occur on faults or segments that have experienced positive CST from prior event(s) (Freed, 2005;Mildon et al., 2019;Stein, 1999;Stein et al., 1997;Verdecchia & Carena, 2015, 2016), conversely earthquakes are less likely to occur in negative 10.1029/2023JB026627 9 of 16 CST regions (stress shadows, Harris & Simpson, 1998).In addition, earthquakes tend to rupture when the CST exceeds a triggering threshold in the range of 0.1-1 MPa (Scholz, 2010).Therefore, we assume that earthquakes should predominantly occur on faults with predominantly positive cumulative CST above a defined triggering threshold.Large earthquakes impact the regional background seismicity, triggering or retarding smaller earthquakes (Scholz, 2010;Toda et al., 2012) and aftershock distributions broadly follow the patterns of positive CST (e.g., King et al., 1994;Toda et al., 1998).Therefore, we assume that for the more plausible earthquake sequences, there will be a higher number of instrumental earthquakes that occur in positively stressed regions.
Based on the arguments above, we compare and filter all modeled sequences by four criteria; We use the median of all models as an initial filter, removing the lower 50% of models and keeping the other half which meet our four criteria more closely.We apply the criteria independently and only keep those earthquake sequence models that exceed the threshold for all criteria.
By evaluating these filtered earthquake sequences, we propose that the rupture scenarios that are physically more plausible will pass all the criteria, have high values, and occur more frequently in the filtered sequences than less plausible rupture scenarios.

Results
Figure 4a shows the results of filtering all 5,184 modeled sequences (each point represents one modeled sequence).
The mean CST on rupture segments prior to rupture is from 0.16 to 0.30 MPa (median 0.22 MPa).The percentage of elements exceeding 0.2 MPa varies between 26.0% and 47.9% (median 35.2%).For elements exceeding 0.5 MPa, the range is 7.1%-16.9%(median 11.5%).The percentage of instrumental earthquakes occurring in positive CST regions ranges between 15.1 and 20.8 (median 17.7%).In general, modeled sequences with lower model numbers tend to lie above the median more often, this is because sequences are combined starting with the preferred scenarios.for example, model 1 features the first scenario of all events, which is the preferred scenario based on available data, while scenarios 2 and 3 are possible, but less likely alternatives.After combining the four filters, 1127 modeled sequences remain (red dots in Figure 4a).For each earthquake with two or three possible rupture scenarios, we calculated the proportions of the input rupture scenarios in the filtered modeled sequences (Figure 4b).For the 1,653, 1,717, 1,920, and 1,928 earthquakes, the input scenarios occur in similar proportions in the filtered modeled sequences.However, scenarios 1 and 2 are very similar for the 1,653 earthquake, as both rupture the same segment with slightly different parameters, and when combined, they outnumber scenario 3 which ruptures a different segment.For some events, one or two scenarios dominate in the filtered models, for example, the 1,645, 1,651, 1,653, 1,880, 1,895, and 1,899 events.The 1,899 event shows a particularly strong indication toward scenario 1.
Figure 4c shows the mean cumulative CST for each fault element from all 1127 filtered sequences in 2020.This gives us a broad idea of the current state of stress on the fault network, given the uncertainties regarding some of the historical earthquake scenarios.The result shows notable positive CST build-up in two areas of the BMGF and on the Pamukkale fault.The western part of the BMGF and the Kuşadası fault also have slightly elevated CST.

Interpretation of Historical Source Faults
Following the approach and assumptions outlined above, we propose a novel way to use Coulomb stress modeling to help constrain the rupture extent of historical earthquakes.
Owing to the method of combining the earthquake scenarios, patterns can be observed in the remaining modeled sequences.For example, the first (1,645) earthquake in the sequence is featured in third of the modeled sequences with each of its scenarios (the first 1,728 sequences feature scenario 1, the next feature scenario 2 etc.).After 1728 modeled sequences, a notable increase in models above the threshold (red) is observable in Figure 4a.This agrees with the dominant occurrence of scenario 2 of the 1,645 earthquake in the remaining models (Figure 4b).In general, earlier combinations (lower numbers in Figure 4a) are more plausible according to the historic record since they feature more primary scenarios, while higher model numbers feature more of the less plausible alternatives.For the 1,899 earthquake, our modeling results favor the first scenario, rupturing the Atça and Pamukören segments (Figure 1c), as this scenario appears in a high proportion of the filtered solutions.A third possible scenario consisting of two separate earthquakes (Ocakoğlu et al., 2013) was not assessed here, but this scenario would be similar to the second scenario, which we determine to be less plausible.For the 1,645 and 1,651 earthquakes, the first two scenarios appear to be more likely than the third; however, as these are the first two earthquakes modeled, our results and interpretation is likely to be biased by the initial state of stress in the model.For most events, the first scenario dominates the filtered models, showing that historical data generally provides good constraints.One exception to this is the 1880 earthquake where scenario 2 (no rupture modeled on the KMGF) dominates, however we do not consider this to be a reliable or conclusive indication that the 1880 earthquake was mis-recorded.This is because the KMGF is across-strike of the BMGF, and therefore it is consistently negatively stressed by earthquakes on the BMGF.For earthquakes featuring equal proportions of all scenarios (especially 1717, 1920, and 1928) it is challenging to say which scenario is most likely, hence perhaps these earthquakes merit further study.However, given that scenario 1 and 2 for the 1,653 earthquake are spatially very similar and together they dominate in the filtered sequences, we suggest that rupture at the eastern end of the BMGF is most likely.Furthermore, this example shows that our approach is less effective when there are minor spatial differences between scenarios.
Considering the model time span (∼400 years), the slip rates of the faults (up to 1.4 mm/a), and the coseismic slip (max.0.8-4 m), recurrence of earthquakes on the same fault segments within this sequence is generally unlikely.This also plays a role in filtering the modeled sequences, reducing the likelihood of models featuring two events on the same segment.Accordingly, if the 1,899 earthquake ruptured the Atça segment (as our models suggest), it is plausible that the 1,645 earthquake occurred on a different segment, probably the Umurlu segment (scenario 2), which also occurs more often in the remaining models (Figure 4b).
Our models exclude the M W 5.6 earthquake occurring in the eastern BMG in 1986, which could be promoted by CST.Epicenter locations place it at the Pamukören segment (Leptokaropoulos et al., 2013;Taymaz, 1993), which possibly experienced a static stress increase due to the 1920 rupture on the adjacent Atça segment, although our results are unclear about its source fault segment.Due to the long time span between both events a stress-triggering relationship is questionable, though an earthquake on the eastern Atça segment would indeed favor rupture of the Pamukören segment.
Figure 4c displays the average cumulative CST of all filtered models in 2020, giving insights into the current stress state and therefore potentially the seismic hazard, with the following caveats.Narrow width zones of increased stress, such as those on the BMGF, are probably a result of how the rupture lengths and exact locations of the ruptured segments were chosen.The assumption that stress is fully released on ruptured segments is uncertain and stress changes from postseismic relaxation are not considered.Therefore, our results can be used to infer general trends, but no detailed seismic hazard interpretations should be drawn.Figure 4c shows a heterogeneous stress state in the Denizli Basin, featuring some highly stressed fault segments.The eastern BMGF and the Kuşadası fault also feature elevated CST.These observations complement patterns of recent seismicity, which is concentrated on the western and eastern ends of the BMGF (Figure 3c).Earthquakes on the BMGF cause negative CST on the KMGF, which may explain the relative seismic quiescence during recorded history and the debated fault activity.We suggest that the KMGF should be considered as active, given the geomorphology, our field observations (Figure 2) and our hypothesis that the fault is in a stress shadow from earthquake activity on the BMGF.
There is debate about the CST triggering threshold, so we tested the sensitivity of changes to the filters and filtering thresholds on our interpretations of the most plausible rupture scenarios (Figure 5).In addition to the four applied filters, we tested filtering with the percentage of elements exceeding 0.0 MPa (Ziv & Rubin, 2000) and 0.1 MPa (Kilb et al., 2002).As with the other filters, all modeled sequences are sorted by the percentage of elements exceeding the threshold and the models below the median are removed.Variations between individually applied filters (Figure 5a) are in a range that the scenario with the highest proportion is almost always unchanged.Special attention should be drawn to the filter based on seismicity in stress-triggering zones (Figure 5a, plot F).
Although it is based on a slightly different modeling approach than the other filters, the scenario proportions generally complement the other results.The overall low percentage of earthquakes in stress-triggering zones are due to earthquake clusters within the assessed map frame but outside the areas influenced by CST.These clusters affect all models to the same extent and hence do not influence the result.Filtering by seismicity in stress-triggering zones is mainly influenced by recent earthquakes, however earlier earthquakes have an impact on these, so the earlier events have an indirect but important impact on this filter.When combining different filters (Figure 5b), the results yield the same trends with minor changes to proportions.We also tested whether the applied filtering threshold of 50% (median) yields different results than higher or lower thresholds (Figure 5c).When changing the quantile of removed models to 33% or 66%, the number of remaining filtered models changes to 2,122 or 457 respectively, but our interpretations of most plausible earthquake scenarios remain similar.
Furthermore, several model parameters can cause minor impacts on the filtering process.Firstly, the used fault geometries are averaged from multiple sources/measurements, usually at outcrops.Variations in fault strike or dip can cause minor changes to the Coulomb stress transferred, though the overall observed patterns of positive and negative stress remain similar (e.g., Mildon et al., 2016).Secondly, the slip rates used to calculate interseismic loading, which are derived from different field-based methods, could contain errors.Absolute slip rates may vary, and the fastest slip may not always be in the fault center.For example, Özpolat et al. (2022) suggest that rock uplift rates along the BMGF are larger in the easternmost segment, whereas we model the highest slip rates in the center of the faults.The actual slip rates for the studied period are difficult to estimate.If the magnitude of the slip rates used for interseismic loading were underestimates, the modeled CST would correspondingly be lower, and thus uncertainties in the slip rates affect the magnitude of the interseismic CST, but not the broad pattern of positive and negative stress.Furthermore, the cumulative CST is dominated by the larger-magnitude coseismic CST.
Thirdly, CST from large earthquakes outside the study area can influence the stress state of the faults and the distribution of microseismicity.Owing to the distance to other faults we assume that these influences are small, however minor variations should not be ruled out.

Applicability to Other Historical Sequences and Tectonically Active Areas
In contrast to previous studies, this new approach presented herein does not simply assess the difference in CST induced by uncertain source fault parameters, but also gives an estimate of the most likely source faults for historical earthquakes, thus improving the reliability of future models and other studies relying on historical earthquake records.This approach could be applied to any historical earthquake sequence that is lacking constraints on some source faults, irrespective of the geological setting or specific earthquake sequence.Although we suggest that the results will be clearer and more useful if at least some of the events in a studied earthquake sequence are already well-constrained from other studies (e.g., paleoseismology, mapped surface ruptures, historical damage patterns), since this reduces the number of possible scenario combinations and gives better constraints in the filtering process.
If this approach was applied to other regions, we suggest that the filters or filter combinations used will probably depend on the region, timescales and magnitude of earthquakes.Using the percentage of seismicity in stress-triggering zones is probably more reliable for earthquake sequences with more recent events, which are likely to have a bigger impact on the recent stress field.We also add caution for utilizing the 0 MPa threshold filter because most elements are positive due to interseismic loading but could be effectively used to remove ruptures on segments experiencing negative CST.

Conclusions
We present a novel use of Coulomb stress modeling to add constraints to historical earthquake sequences.We introduce a new version of the 3D-faults code, originally published by Mildon et al. (2016).This is used to model the Coulomb stress (coseismic and interseismic) evolution over 400 years for the Büyük-Menderes Graben region.By modeling multiple possible historical earthquake scenarios and subsequent evaluation of >5,000 earthquake sequences, we contribute to the search for the source fault for several events in the region, including the debated 1899 event.We propose that the 1899 event is mostly likely to have occurred on the Atça and Pamukören segments.Combining the most plausible scenarios, we have calculated the average state of stress on all faults in the region, which implies that the Pamukkale fault is highly stressed across the range of models we test, in agreement with the high seismicity in the area.
In general, the clarity and robustness of results for each earthquake depends on the available data for all events in the sequence.If source faults for most events are well constrained, the method is more likely to give a robust result for a single, poorly constrained earthquake, whereas sequences with sparse knowledge on most events will yield less certain results.Additional evaluation for other earthquake sequences would be beneficial to test the method's capabilities and limitations.significantly improved the quality of the manuscript.Thanks to Alessandro Verdecchia for guidance and discussion on using the backslip method.M.D. is supported by a University of Plymouth School of Geography, Earth, and Environmental Sciences PhD studentship.Original development of the 3D-faults code was supported by NERC PhD Studentship NE/L501700/1 and JSPS Short Term Fellowship PE15776 to ZM.For questions and comments concerning the '3D-Faults' code please contact M.D. or Z.M. E.H. publishes with permission from the Executive Director of the British Geological Survey.
evaluate possible source faults of destructive historical earthquakes using Coulomb stress modeling • Results constrain the extent of historical ruptures in W Türkiye, particularly the 1899 earthquake which affected a ∼100 km wide area • Areas of positive Coulomb stress in the Büyük Menderes Graben and Denizli Basin broadly correlate with instrumental seismicity Supporting Information: Supporting Information may be found in the online version of this article.

Figure 2 .
Figure 2. Field data constraining some of the fault geometries used for the model (Küçük Menderes Graben fault and Honaz fault).(a) Location of the study area in western Türkiye.(b) Map of outcrop locations (red dots).Panels (c)-(e) show outcrop photos and structural data as lower hemisphere projections of fault planes (great circles), striations (black dots) and lineaments deduced from corrugations (red diamonds).Mean values are given as dip direction/dip (planes) and trend/plunge (lineaments).(c) Outcrop of the Küçük Menderes Graben fault featuring two small but well-exposed and fresh fault planes with striations indicating normal faulting.(d) Fault plane of the Honaz fault near Honaz town center.Note that only 3 striations could be measured, hence lineaments from corrugations are more reliable indicators of slip sense.(e) Large, corrugated fault plane of the Honaz fault, ∼1 km west of Honaz town center.
Pamukören (S2) segment, although other scenarios are possible.In 1928 another poorly constrained earthquake was recorded in the KMG (KOERI, 2020).Similar to the 1880 event we model one scenario on the KMGF and a second leaving the event out.The 1955 Söke-Balat earthquake (N.N.Ambraseys & Jackson, 1998;Eyidogan & Jackson, 1985;Gürer et al., 2009;Paton, 1992) ruptured the Söke (Priene-Sazlı) fault (one scenario).A M W 5.6 earthquake likely ruptured the Pamukkale fault in 1963(KOERI, 2020).We only model one scenario in the central segment but the actual segment that ruptured is unknown.The last earthquake in our model is the 1965 M W 5.7 earthquake (KOERI, 2020) which, according toKumsar et al. (2016), ruptured the Honaz fault (one scenario).

Figure 3 .
Figure 3. Workflow of combining scenarios and modeling Coulomb Stress Transfer (CST) throughout earthquake sequences.(a) Using a logic tree approach to combine 1-3 input scenarios for the 16 modeled earthquakes, yielding 5,184 possible earthquake sequences.(b) An example of cumulative (interseismic plus coseismic) CST resolved on the fault network after the 1899 earthquake.(c) Cumulative CST resolved on specified receiver faults 180°/60°/−90° (dip direction/dip/rake, depth 8 km) for one earthquake sequence.Black dots indicate recorded seismicity (ISC, 2022) after the last event (1965-2020).
(a) mean cumulative CST on the source fault prior to rupture, (b) percentage of fault elements exceeding 0.2 MPa and (c) 0.5 MPa, and (d) percentage of earthquakes (M W > 3.5 from 1965 to 2020 (ISC, 2022)) occurring in positively stressed areas.

Figure 4 .
Figure 4. Evaluating the modeled earthquake sequences.(a) Filtering of modeled sequences.Each dot represents one earthquake sequence.(b) Proportions of input scenarios in the models remaining after filtering.(c) Average Coulomb Stress Transfer on the fault network from all 1,127 filtered models.

Figure 5 .
Figure 5. Testing the sensitivity of the resulting scenario proportions on different filtering.(a) Independent application of six different filtering metrics.While the filters A-E are based on Coulomb Stress Transfer (CST) resolved onto the fault network, filter F uses the CST resolved on specified receiver planes, which is then compared to background seismicity on a map.(b) Combining filtering metrics.A combination of filters A, D, E and F is used in our study (Figure4).The depicted alternatives show minor variation and do not change which scenario is the frequent.(c) Using the same filter combination as in the study but varying the filtering threshold.In the study the median was used (center plot).Scenario proportions and the resulting interpretations regarding source faults change negligibly for a 33%-threshold but more noticeably for the 66%-threshold (removing 2/3 of models by each filter).

Table 2
Model Parameters of All Earthquake Scenarios