Under-Displaced Normal Faults: Strain Accommodation Along an Early-Stage Rift-Bounding Fault in the Southern Malawi Rift

One of the fundamental problems in continental rift segmentation and propagation is how strain is accommodated along large rift-bounding faults (border faults) since the segmentation of propagating border faults control the expression of rift zones, syn-rift depo-centers, and long-term basin evolution. In the Southern Malawi Rift, where previous studies on the early-stage rifting only assessed border fault structure from surficial and topographic expression, we integrate surface and subsurface data to investigate border fault segmentation, linkage, and growth as proxies for strain accommodation along the Bilila-Mtakataka Fault (BMF) System. We used 30 m-resolution topographic relief maps, electrical resistivity tomography (ERT), and high-resolution aeromagnetic data to characterize the detailed fault geometry and provide a more robust estimate of along-fault displacement distribution. Our results reveal a discrepancy between sub-aerial segmentation of the BMF geometry (six segments), scarp height (five segments) reflecting the most recent episodes of fault offset, and cumulative throw (three composite segments) reflecting the long-term fault offset. We also observe that although the BMF exhibits continuity of sub-aerial scarps along its length, the throw distribution shows a higher estimate at the Northern-to-Central segment relay zone (423 m absolute, 364 m moving median) compared to the Central-to-Southern segment relay zone (371 m absolute, 297 m moving median). The ERT profiles across the relay zones suggest a shallower basement and a possible canyon-mouth alluvial fan stratigraphy at the Central-to-Southern segment relay zone, contrasting the deeper basement and “simpler” electrical stratigraphy at the Northern-to-Central relay. The results suggest a more complex long-term evolution of the BMF than was assumed in previous studies. A comparison of BMF’s maximum displacement-vs-length with those of other Malawi Rift border faults and global normal fault populations suggest that although the BMF has possibly reached its maximum length, it remains largely under-displaced as its 580–837 m maximum displacement is significantly lower than that of faults of equivalent length. We suggest that the BMF may continue to accrue significant strain as tectonic extension progresses in the Southern Malawi Rift, thus posing a major seismic hazard in the region.


INTRODUCTION
Normal faults grow as tectonic extension causes them to slip. During continental rifting, normal faults grow by strain accommodation leading to the progressive evolution of their associated rift basins Nixon et al., 2016;Muirhead et al., 2019). Strain accommodation, strain transfer, and deformational styles in early-stage continental rifts are often influenced by regional tectonics, stress orientation, lateral variation in crustal stretching rates, pre-existing basement structures, lithology, magmatism, and lithosphere rheology (e.g., Nixon et al., 2016;Peacock et al., 2017;Muirhead et al., 2019;Rotevatn et al., 2019;Wright et al., 2020;Kolawole et al., 2021a;Gouiza and Naliboff, 2021). During continental extension, strain is accommodated by crustal deformation in the form of rift faulting, and the distribution and pattern of faulting can be quantified from a series of structural and tectonic attributes associated with fault growth (Bonini et al., 2005). These attributes include the younging direction of the border faults (Abbate and Sagri, 1980;WoldeGabriel et al., 1990;, the changes in border fault length and the elastic thickness of the lithosphere (Hayward and Ebinger, 1996), lateral variation in ages of syn-rift deposits (Maguire et al., 2003;Furman et al., 2004), and lateral migration of magmatic activity (Zanettin et al., 1979).
In rift basins, upper crustal extension is accommodated by two groups of fault systems: border faults and intra-rift faults (e.g., Corti, 2009). The border faults are the largest faults bounding a rift segment and often accommodate the most offset, whereas intra-rift faults are the smaller faults that are found within the inner depressions of rift basins (Scholz and Contreras, 1998;Goldsworthy and Jackson, 2001;Ebinger, 2005;Corti, 2009;Muirhead et al., 2019). However, the new and recent 2015 multidisciplinary Study of Extension and maGmatism in Malawi and Tanzania (SEGMeNT) project provided seismic data data over Lake Malawi that show some very long intra rift faults with large throws of several km, even larger than the throw of some border faults in Southern Malawi Rift Shillington et al., 2016;.
Large normal faults typically develop by an initial nucleation of isolated segments, which progressively propagate and link up with one another. The interaction and linkage of adjacent propagating fault segments typically create structurally complex areas, referred to as "strain transfer zones" or "relay zones" (e.g., Morley et al., 1990;Morley, 1999;Jackson and Rotevatn, 2013;Childs et al., 2017). These relay zones accommodate the along-and across-axis variations in the magnitude of subsidence of grabens and the elevation of uplifted rift flanks or normal fault footwall blocks (Rosendahl, 1987;Morley et al., 1990;Faulds and Varga, 1998;Morley, 1999).
There exists an empirical relationship between fault length and displacement, which provides insight into the mechanics of fault growth. It presents a useful tool for characterizing the extent and rates of tectonic deformation in continental rift zones (Bergen and Shaw, 2010). Scaling laws of maximum displacement versus fault length exhibit a power-law relationship (e.g., Walsh and Watterson, 1988;Cowie and Scholz, 1992a;Cowie and Scholz, 1992b;Dawers et al., 1993;Manighetti et al., 2001). Normal fault propagation models suggest that faults grow along fault tips, and that there is a proportional increase in the length of the fault and maximum displacement throughout the slip history of the fault ("propagating fault model", e.g., Walsh and Watterson, 1988;Dawers et al., 1993;Cartwright et al., 1995;Walsh et al., 2003). Another model ("constant-length model") suggest that the total length of the fault is attained early in the fault history and that most displacement is accumulated after the total length has been established (Walsh et al., 2003;Nicol et al., 2005;Densmore et al., 2007;Schlagenhauf et al., 2008;Jackson and Rotevatn, 2013;Nixon et al., 2016). Rotevatn et al. (2019) proposed a hybrid fault growth model which suggests that normal faults initially grow by lateral propagation in the early stages of the fault lifespan (20-30% of the fault lifespan) to establish the total length of the fault, then, subsequently grow by constant-length model for the rest of the fault lifespan (70-80% of the fault lifespan). Also, recent numerical models by Pan et al. (2020) revealed that normal fault arrays evolve via alternating phases of faultlengthening and localization. However, the process of strain accommodation by normal fault growth during the early stages of continental extension remains incompletely understood. Since the segments of important rift faults are often established during early-stage rifting, this knowledge-gap limits the complete understanding of the processes that govern the variation of structural architecture along margins of continental break-up.
In this contribution, we investigate the detailed structure and distribution of displacement along the Bilila-Mtakataka Fault (BMF), an active border fault in the Southern Malawi Rift, East Africa ( Figures 1A-C), to understand the pattern of fault growth and strain accommodation in areas of early-stage continental extension. This fault is suggested to be seismogenic and may have hosted the M 6.2 earthquake, 9 km NNW of Salima, Malawi at a depth of about 30.3 km (from USGS earthquake hazard program). In addition, the results of the study provide insight into the seismic hazard of large riftbounding normal faults.

Geologic and Tectonic Setting of the Malawi Rift
Eastern Africa has witnessed three phases of tectonic extension during the Phanerozoic. The East African Rift System (EARS) represents the third and currently active phase of extension which began in the Oligocene (Delvaux, 1991;Ebinger and Scholz, 2011;Roberts et al., 2012;Ebinger et al., 2019). Thẽ 750 km-long N-S trending Malawi rift represents a prominent segment of the EARS (Ebinger et al., 1987;Ebinger et al., 1989), and is located near the southern tip of the magma-poor western branch of the rift system ( Figures  1A,B). The Malawi Rift is sub-divided into eight distinct 100-150 km-long grabens and half grabens which exhibit alternating polarities along-strike (Laó-Dávila et al., 2015). Using a complete dataset that incorporated seismic imaging across Lake Malawi and estimation of fault throws and displacements, Scholz et al. (2020) provide an alternative and simpler structural framework than Laó-Dávila et al. (2015) for the Malawi Rift where only three main segments are presented. To the north, the northernmost graben, the "Karonga" or North Basin terminates at the Rungwe Volcanic Province, whereas, the southernmost graben of the Malawi Rift is defined by the Zomba Graben which terminates at a zone of exposed basement extending into the Shire Rift Zone (Kolawole et al., 2021b).
Although the Malawi Rift largely developed during the Cenozoic evolution of the EARS, it is proposed to have propagated southwards based on: 1) the extensive outcrops of Cenozoic syn-rift sequences in the northern basins of the rift which are absent in the south (Delvaux, 2001), and southward thinning of the syn-rift sedimentary cover (Specht and Rosendahl, 1989;Flannery and Rosendahl, 1990), and 2) southward decrease in the throw of border faults with corresponding changes in the seismic depositional facies of the syn-rift sequences . However, it is also plausible that the initiation of Cenozoic rifting across the entire length of the Malawi Rift could have been synchronous and that rifting occurs at different rates along the rift length (Ojo et al., 2021). Overall, crustal stretching rates decrease from north (4.5 mm yr −1 ) to south 1.3 mm yr −1 . The major graben and some associated major border faults in the Southern Malawi Rift are the Makanjira Graben (Bilila-Mtakataka and Malombe faults), and the Zomba Graben (Zomba and Lisungwe faults).

The Bilila-Mtakataka Fault
The Bilila-Mtakataka Fault (BMF) is one of the most prominent rift-bounding faults in Southern Malawi Rift, characterized by a steep, clean fault scarp of about 4-15 m height, and was first highlighted in Walshaw (1965) and Dawson and Kirkpatrick (1968)  Although the fault has been estimated to be an outstanding >100 km-long continuous normal fault scarp (Jackson and Blenkinsop, 1997;Hodge et al., 2018), the subsurface structure of the fault remains unknown. Based on the large length of the BMF, Jackson and Blenkinsop (1997) proposed the possible control of a thick seismogenic layer on the fault evolution. The BMF is suggested to have developed within a strong lithosphere with an effective elastic thickness of about 30 km (Jackson and Blenkinsop, 1993).
The surficial geomorphology of the BMF comprises a metamorphic basement that is composed of fractured and weathered basement rock at near-surface but still intact at depth. It is also made up of a deeper fresh basement, above which the syn-rift sediments accumulate in the hanging wall of the BMF. The sediments are primarily alluvial sediments with colluvial lobes confined to the foot of the fault escarpment (Walshaw, 1965). The sedimentary basin in the hanging wall is a fluvio-lacustrine depositional environment with occasional sediment inundation from 100-year cycles of over-flooding events in Lake Malawi (Dulanya et al., 2014;Dulanya, 2017).
The fluvio-lacustrine depositional environment is due to several transverse cutting streams that originate from the elevated footwall areas (e.g., Nsipe-Livelezi Shelf watershed) and flow downslope, across the BMF escarpment into the basin where they either empty into Lake Malawi or join the Bwanje River which is the prominent axial stream draining the Makanjira Subbasin ( Figure 2D). The northward-flowing Bwanje River originates from a watershed in an elevated area near Balaka that marks the southern tip of the BMF, separating the Bwanje valley in the Makajinra Trough from the Liwawadzi Valley in Zomba Graben to the south (Walshaw, 1965). Hodge et al. (2018) investigated the sub-aerial map-view geometry of the BMF trace and controls of basement fabrics using high-resolution satellite topography data and field measurements. Based on the map-view along-trend changes in fault geometry and the along-fault distribution of vertical surface separation, six segments were identified. These segments consist of Ngodzi, Mtakataka, Mua, Kasinje, Chitsulo, and Bilila geometrical segments. The BMF was suggested to have a scarp height distribution that is likely influenced by the map-view fault geometry, and that the fault geometry developed as an exploitation of the local metamorphic foliations (Hodge et al., 2018;. In this contribution, we describe the surface expression of the subaerial BMF trace as a representation of its map-view geometrical segmentation, using the names given by Hodge et al. (2018). In this study, we integrate topographic analysis with geophysical imaging of the subsurface structure along the Bilila-Mtakataka Fault (BMF). We obtain along-fault scarp height and footwall height measurements from one arc-second (30 m spatial resolution) Shuttle Radar Topography Mission (SRTM) elevation data. Scarp height provides the sub-aerial component of recent fault throw, and footwall height provides a measure of sub-aerial component of long-term cumulative fault throw. We utilize aeromagnetic transform techniques to generate a depth-tomagnetic basement map which provides estimates of the subsurface component of long-term cumulative fault throw. Our sampling points from both the SRTM topographic data and aeromagnetic depth-to-basement are approximately the same along the entire trace length of the fault. At each sample point, the estimated throw represents its relative value across a fault surface, i.e., relative to depth, or its maximum value (Tao and Alves, 2019). Profiles were taken perpendicularly across the fault at every 1 km because the sampling interval/fault length ratio (δ) has to be <0.05 for adequate resolution and to avoid loss of vital information. The δ = 0.05 threshold value is recommended for faults of any scale and profiles should be taken at intervals that are <5% of the total length of fault if they are believed to be segmented (Tao and Alves, 2019). We utilize Electrical Resistivity Tomography (ERT) to investigate the fault zone and assess the lateral variation of local hanging wall stratigraphic architecture along the trend of the fault.

Scarp and Footwall Height Measurements From Digital Elevation Models
The fault scarp indicates the most recent slip along a normal fault. The scarp heights were extracted and measured from a 30-m resolution Shuttle Radar Topography Mission (SRTM)-Digital Elevation Model (DEM). Fault scarps were outlined and traced from a hillshade map created on QGIS using the SRTM-DEM. The scarp height (h) is the vertical distance of the steepest slope along the local footwall elevation. The vertical distance of the local footwall slope is the local footwall elevation, i.e., the top of the local footwall slope to the top of the local hanging wall within which the fault scarp is hosted (h'; Caskey, 1995;Hornsby et al., 2020). We measured the vertical distance of the footwall between the steepest points along the fault scarps (as shown in Figure 3). The scarp heights and local footwall heights were measured and plotted against fault length to produce a displacement-length relationship. The methodology adopted in this study accounts for variables such as ancient topography, weathering and erosion, and structural irregularities due to the subjective nature of this process and observed irregularities along the fault trace. Our measurements are likely to be underestimated as the escarpment continues up to a higher elevation, however, since it is difficult to constrain the component of paleotopography in the footwall elevation, we take a conservative approach.
We note that the method of scarp height measurement adopted in this study is different from the methodology used by Hodge et al. (2018; which adopted a semi-automated approach for extraction of relief height from satellite topography data. Their study measured the "vertical separation" as an estimate of fault throw rather than scarp height. Vertical separation (or surface offset) is distinguished from scarp height in literature (Bucknam and Anderson, 1979). Further, to better capture the full length of the BMF and understand the long-term fault growth, we extended our study to include the northern extents of the BMF fault trace than the extent covered in Hodge et al. (2018). Segments on the scarp height plot along the fault length were classified as part of the fault between two consecutive minima separated by the maximum point between them. Relay zones on the along-fault scarp height plot are classified as flat parts of the fault between two consecutive minima and not separated by any maximum points.

Estimation of Depth-to-Magnetic Basement From Aeromagnetic Data
To assess the along-fault distribution of subsurface components of fault throw, we estimated the depth to the top of the magnetic basement (d) in the hanging wall of the BMF using 62 m spatial resolution aeromagnetic data from Malawi. The aeromagnetic data were acquired in 2013 (source: Geological Survey Department of Malawi) with 80 m terrane clearance along NE-SW lines with a line spacing of 250 m. Since the primary and dominant magnetic source in the study area is the Precmabrian gneissic basement ( Figure 2B), d is taken to be a measure of the depth to the crystalline basement. For our basement depth calculation, we use the Source Parameter Imaging (SPI) transform of the total magnetic intensity (TMI) anomaly grid (Blakely and Simpson, 1986;Thurston and Smith, 1997;Smith et al., 1998).
The SPI formulation assumes a step-type source in which the depth to the top of the magnetic source is an inverse of the peak value of the local wavenumber over the step source.
where K max is the peak value of the local wavenumber K over the step source, T is the Tilt derivative, VDR is the first vertical derivative, THDR is the total horizontal derivative, and M is the total magnetic field anomaly (TMI). The computation of the SPI depth solutions across a TMI grid involves an initial computation of the local wavenumber from the tilt derivative, then smoothing of the wavenumber grid using a Hanning filter to reduce noise, and finally, localized peak wavenumber detection using a Blakely test (Blakely and Simpson, 1986). Finally, we extracted the values of depth-tomagnetic basement in the hanging wall at a distance of 1 km basinward from the base of the fault scarp at each scarp height Frontiers in Earth Science | www.frontiersin.org April 2022 | Volume 10 | Article 846389 sample point ( Figure 3). We use the 1 km extraction distance to maintain consistency, since the deepest sections near the BMF segment centers, and shallowest sections in the relay zones commonly occur within this distance from the base-scarp. The estimation of depth-to-magnetic basement from aeromagnetic data has an inferred accuracy of ±20% (Gay, 2009). Thus, we indicate an error bar of ±20% of the measurements in the depthto-basement scatter plot.

Estimation of Cumulative Fault Throw and Maximum Displacement along the Bilila-Mtakataka Fault
Fault throw is the vertical component of fault offset between correlative footwall and hanging-wall of the fault (Muirhead et al., 2016;Tao and Alves, 2019). Thus, we estimate the cumulative fault throw as the sum total of both the surface and subsurface components of throw along the fault. To obtain the minimum throw (T) at each sample point along the BMF, we use: where h′ is the measured local footwall height, and d is the depthto-magnetic basement. Along-fault throw-distance plots present a view of fault segmentation which then provides insights into the complex structural evolution of the BMF in terms of its linkage and long-term growth. Therefore, the measured throws along the BMF serve as proxies for the overall displacement along the fault. However, to calculate maximum dip-displacement along the BMF fault plane D, we use dip magnitudes θ measured directly from the electrical resistivity tomography images (section 3.4) and the estimated maximum throw T max following the relation: Based on pure normal faulting focal mechanism solutions along the BMF (Stevens et al., 2021), we assume a faultorthogonal slip vector for the BMF. In addition, the lateral projection of clusters of microseismicity at depth beneath the BMF delineates a dip of~42° (Stevens et al., 2021), thus expanding the range of our estimated maximum displacement for the fault, considering the possibility of non-planarity of faults at depth. To better understand the long-term growth of the BMF, we compare its displacement-length ratio with those of other well-studied Malawi Rift border faults (Livingstone, Usisya, and South Basin border faults; Figure 1B). For these other border faults, we obtain estimates of their dip, sub-aerial components of throw and dipdisplacement from previously published seismic reflection images along Lake Malawi, assuming a normal slip vector for the faults (Specht and Rosendahl, 1989;Scholz et al., 2020;Shillington et al., 2020). The cumulative throws and displacements for these other border faults were estimated by adding the sub-aerial components from seismic data to the surface components of the fault escarpments.

Electrical Resistivity Tomography of Fault Zone and Local Hanging Wall Stratigraphic Architecture
The BMF is composed of multiple segments as defined by surface scarp height variation and plan-view geometry (this study and Hodge et al., 2018). To characterize the local subsurface electrical structure in the hanging wall of the BMF, we acquired four 2-D electrical resistivity profiles across the segments and their relay zones ( Figure 2A). Two of the survey profiles were acquired near the central part of two of the fault segments (ERT Profiles 1 and 4), and the other two profiles were acquired near two of the relay zones (ERT Profiles 2 and 3). The survey transects are oriented perpendicular to the fault trend, extending across the BMF fault scarp at each of the survey locations, except for a section along ERT Profile 3 which is slightly oblique to the fault trend due to land accessibility challenges in the field ( Figures 2C-E). Along ERT Profile 3, the fault-oblique section occurs across the fault scarp, whereas the southwestern and northeastern sections of the profile are at higher angles to the fault trend ( Figure 2D). At all the survey locations, the profiles extend from the footwall of the fault, across the top and base scarps, and terminate at a considerable distance into the hanging wall ( Figures 2F,G). We utilized a 10 channel Iris Syscal Pro resistivity meter with 72 electrodes and a 10 m electrode spacing along the four profiles ( Figures 2H,I). Due to field constraints, three of the profiles (ERT Profiles 1, 2, and 4) were each 710 m-long, and the fourth (ERT Profile 3) was 1,070 m-long (acquired using roll-along sequence). Based on the ERT profile lengths, we consider that the surveys only image the fault zone and electrical stratigraphy of the local hanging wall slope of the fault. Except for ERT Profile 1 at which only the dipole-dipole array was acquired, we acquired both Wenner-Schlumberger and dipole-dipole arrays for all the profiles. The Wenner-Schlumberger array is sensitive to both horizontal and vertical structures, while the dipole-dipole array is good for mapping vertical structures (Loke, 2002).
To obtain a dense dataset with good lateral and vertical resolution, we combined the data from both array types to create a mixed array using Prosys II (Zhou, et al., 2002). Bad data points were culled using the RMS statistics based on preliminary inversion, and the final inversion model results have RMS error <10%. Supplementary Table S1 contains the relevant details of the inversion parameters applied. The data processing and inversion were performed using RES2DINV software (Loke and Barker, 1996) using the smoothnessconstrained least-squares method (Claerbout and Muir 1973;deGroot-Hedlin and Constable 1990). The least squares formulation provides a RMS (root-mean-squared) value which is the factor of discrepancy that describes the differences between the logarithms of the measured and calculated apparent resistivity values, such that the RMS error value of a model inversion is the measure of data fit for the inversion. The parameterization of the inversion carried out in our study is provided in Supplementary Table S1 of the supplementary document. We present the results of the 2-D resistivity tomography along with plots of model sensitivity estimates.
The model sensitivity is a measure of the reliability of the inversion model (Loke, 2002). The higher the sensitivity value, the more reliable the model resistivity value and vice-versa (Loke, 2002). We defined low (less than 0.1), medium (values between 0.1 and 2.0), and high (greater than 2.0) sensitivity values for each inversion model (Supplementary Figures S1, S2) and used the model sensitivity to constrain the interpretation of the resistivity sections. Low-sensitivity regions are the shaded polygons in the ERT inverse models (Figures 7-10), and are given less weight in the interpretations.
Although there is no borehole data available in the study area to provide a strong constraint on the inversion profile interpretations, we adopt an "in-context" approach for interpretation, which relies on the tectono-geomorphic setting of the target structure (NE-dipping active normal fault with exposed basement footwall uplift; Figure 2A), the surface locations of the exposed crystalline basement and sedimentary cover along the survey profiles ( Figure 2B), and the depositional environment of the basin (humid fluvio-lacustrine setting; Dulanya et al., 2014;Dulanya, 2017). Overall, the typical resistivity values for crystalline basement have a range of 1,000-100,000 Ω-m, and those of unconsolidated sedimentary deposits and weathered basement rocks range 2-1,000 Ω-m (Palacky, 1987;Minsley et al., 2010). We describe <10 Ω-m layers as low-resistivity units,~10-600 Ω-m as low to moderate-resistivity units, and~600 to >3,000 Ω-m as moderate-to high-resistivity units.

Along-Fault Distribution of Scarp Height (h) and Local Footwall Height (h')
At the surface, the NW-trending BMF trace is approximately 152 km in length and exhibits plan-view geometrical ( Figure 2A) and scarp segmentation (Figure 4). Based on the along-fault distribution of scarp height and local footwall height, the profiles ( Figure 4) show that the border fault is subdivided into five scarp segments (Scarp Segments 1-5). Scarp Segment 1 is the northernmost segment along the BMF and it is 57 km long and has a maximum scarp height of 31 m and a maximum local footwall height of 61 m, which clearly contrasts the planview geometrical segmentation of this section of the fault (Ngodzi We note that the relay zone (zone of transition between segments) between Scarp Segments 4 and 5 shows the greatest magnitude of scarp height minima compared to the other relay zones. Also, we note that in both plan-view geometry ( Figure 2) and scarp height distribution (Figure 4), the longest "continuous" segments along the BMF are the northern-most segment (Scarp Segment 1) and the southernmost segment (Scarp Segment 5). In between these two long segments, three smaller segments (Scarp Segment 2, 3 and 4) are found.
Overall, segmentation is generally consistent with the subaerial fault trace geometry (from Hodge et al., 2018) except for scarp height segment 1 where our scarp height does not show a pronounced minima between Ngdozi and Mtakataka segments. However, our maximum footwall height distribution shows a slight minima between Ngodzi and Mtakataka segments and based on these we subdivided scarp height segment 1 into scarp height segments 1a and 1b. The broad trends illustrated by the moving median curves suggest that these five sub-aerial segments can be grouped into three large composite segments considering the relative magnitudes of scarp height minima along the profile. These composite segments consist of a northern composite segment (includes Ngodzi and Mtakataka geometrical segments), a central composite segment (Mua, Kasinje, and Chitsulo geometrical segments), and a southern composite segment (Bilila geometrical segment). The scarp height and local footwall height plots show similar but higher magnitudes in both the northern and central composite segments and appear to gently decrease towards the southern end of the border fault system. Whereas, the maximum local footwall height shows the highest values in the central composite segment and somewhat similar magnitudes in the northern and southern composite segments.

Along-Fault Distribution of Depth-to-Basement
The calculation of depth-to-magnetic basement from the total magnetic intensity aeromagnetic grid ( Figure 5A) using the SPI method produced a grid map of the depth solutions ( Figure 5B). In general, the SPI map shows <500 m basement depths which increases westward towards the BMF escarpment and shallows eastward towards the Shire Horst to the southeast and Lake Malawi to the northeast ( Figure 5B). The shallowing of the basement is most pronounced in the southeast, along the Chitsulo and Bilila segments where the width of sediment cover in the BMF hanging wall is narrow (<7 km-wide) compared to 8-28 km-wide regions further north along the Mtakataka-to-Kasinje segments. We note the presence of a prominent narrow region of shallow basement (~120 m) that extends southwest from the Shire Horst into the relay zone between the Kasinje and Chitsula segments of the BMF ( Figure 5B). In comparison to the other relay zones along the BMF, this relay zone exhibits the shallowest basement. In addition, the narrow region of shallow basement marks the transition from the southern BMF hanging wall region with narrow sedimentary cover to the northern BMF hanging wall region with wider sedimentary cover ( Figure 5B). Along the BMF, the extracted basement depth values show a distribution that roughly mimics the general trend of the maximum local footwall height ( Figure 6). The values are highest near the center but decrease to~180 m near the northern fault tip and~100 m near the southern tip. We observe two prominent minima along the moving median curve: one at distance~42 km, and the other at distance 100 km, thus delineating three large segments which are consistent with the composite segments of the scarp height fitting curve except for the Northern Segment. We refer to these large segments as the Northern Segment (extending between distance 0-40 km; 0-57 km for scarp height fitting curve), Central Segment (distance 45-100 km), and the Southern Segment (distance 105-150 km). The basement attains a maximum depth of~460 m in the Northern Segment,~500 m in the Central Segment, and~460 m in the Southern Segment ( Figure 6).

Along-Fault Distribution of Cumulative Throw, T
We estimated the minimum cumulative throw, T, at each sampling point to produce a distance-throw plot along the BMF by adding the maximum local footwall height and SPI depth-to-basement results ( Figure 6). The fitting curve of the moving median of minimum cumulative throw shows a similar trend as those of both the depthto-basement and maximum local footwall height, consistent with the delineation of three large segments along the BMF ( Figure 6). Overall, the maximum T along the fault is~560 m, located near the center of the entire fault system. The minimum T estimates arẽ 120 m at the southern tip of the fault and~190 m near the northern fault tip. However, the Northern Composite Segment shows a maximum T of~500 m, the Central Composite Segment~560 m, and the Southern Composite Segment~480 m. We consider that the smaller lobes on the throw-distance distribution plot could be data artifacts from either the SPI transform, the moving median fitting curve, or may be possibly representative of early-growth subsegments of the fault, and that the broader lobes are more representative of the long-term fault growth.

Subsurface Electrical Resistivity Structure
In general, the ERT Profiles 1 to 4 (Figures 7-10) commonly show three geoelectrical layers that are truncated and offset by sub-vertical discontinuities. On each profile, the most prominent among the sub-vertical discontinuities is colocated with the base of the BMF escarpment at the surface. The relative sensitivity of the inverse models is mostly composed of model blocks of medium and high sensitivities, particularly within the upper 40 m interval ( Supplementary Information S1). Also, the regions with high resistivity values compared to the surrounding regions generally show low sensitivities (Supplementary Information S1).
ERT Profile 1 was acquired near the center of Scarp Segment 1 (Ngodzi) in the north-western section of the BMF ( Figure 2C) and shows three geoelectric layers (Figure 7). The topmost layer is a~10 m thick laterally-continuous layer dominated by low-to moderate-resistivity values (~10-100 Ωm). This layer overlies a 120 m-thick layer that is dominated by generally low-resistivity zones (<10 Ω-m). The third and basal geoelectric layer appears to be a moderate-to high-resistivity layer (>1,000 Ω-m) that extends up to the near-surface in the northwestern part of all the profiles where it is confined to the footwall of the BMF escarpment. However, we note that the model sensitivity is low in the basal part of the model. ERT Profile 2 ( Figure 8) was acquired near the zone of segmentation (relay zone) between Scarp Segments 1 and 2 (transition from Mtakataka to Mua geometrical segments) in the northwestern section of the BMF ( Figure 2C). The inversion model shows three geoelectrical layers with a 20-40 m-thick topmost layer dominated by~10-300 Ω-m resistivity anomalies, which appear to mimic the electrical character of the topmost layer in ERT Profile line 1. This topmost layer is underlain by a~100 m-thick lowresistivity (<10 Ω-m) layer, which is underlain by a basal layer of elevated resistivity values (>30 Ω-m). Similar to ERT Profile 1, the footwall of the BMF scarp is dominated by high-resistivity anomalies (>1,000 Ω-m).
ERT Profile 3 (Figure 9) was acquired at the relay zone inbetween Scarp Segments 3 and 4 (transition from Kasinje to Citsulo geometrical segments) in the central section of the BMF (Figure 2). The inversion result of this profile contrasts with the other ERT profiles in that it shows a complex geoelectrical model in which the top of the basal high-resistivity layer defines a gently-dipping (not  sub-vertical) slope that is colocated with the base of the BMF scarp. The top-most geoelectric layer is characterized by low-to moderateresistivity (10-60 Ω-m) and has a thickness of 20-40 m which thickens to at least 70 m basinward (near the far eastern end of the profile). This shallow layer overlies a thicker interval (~110 m) defined by discrete bodies of moderately resistive bodies (60-300 Ωm) hosted within a~60 Ω-m background resistivity anomaly. Similar to other profiles, the basal geoelectric layer is composed of high resistivity values (>1,000 Ω-m), extends into and dominates the footwall of the BMF. ERT Profile 4 (Figure 10), acquired near the center of the Scarp Segment 5 (Bilila geometrical segment) in the southernmost section of the BMF, is generally characterized by lowresistivity zones (<30 Ω-m) in the hanging wall of the BMF scarp and very high-resistivity zones in the footwall of the scarp. The BMF scarp is colocated with a strongly subvertical discontinuity that extends across the inversion model.  The inversion results for ERT Profiles 1 and 2 show a 75°dip angle for the BMF electrical gradient (Figures 7, 8), whereas a 70°dip angle is observed on ERT Profile 4 ( Figure 10). These are useful estimates of shallow crustal dip of the BMF, considering that the profiles are orthogonal to the fault trend. Based on a 560 m maximum throw along the BMF, and dips of 70-75°from ERT and 42°from deep microseismicity clusters, we estimate the maximum dip displacement to range between 837 m (upper bound), 580 m (lower bound), and a median of 596 m ( Figure 11). Further north along the Malawi Rift ( Figure 1B), the Livingstone border fault records a maximum displacement of 7.4 km, Usisya border fault records~7.3 km, and South Basin Border Fault records~3.85 km (Figure 11; Specht and Rosendahl, 1989;Scholz et al., 2020;Shillington et al., 2020).

Integration of Topographic and Subsurface Observations Along the Bilila-Mtakataka Fault Zone
The results of this study reveal a distinct difference in the segmentation pattern observable along the surface trace of the BMF (i.e., scarp and geometrical segmentation; Figure 4) and the segmentation pattern shown by the along-fault distribution of depth-to-basement ( Figure 6). Faults usually initiate as isolated fault segments which evolve over time through fault growth and linkage of the isolated fault segments into larger faults (Jackson and Rotevatn, 2013;Childs et al., 2017). As the faults grow through linkage and interactions between the isolated segments, the magnitude of displacement and throw along the resulting fault become cumulative and change over time. The longer the time that has passed since isolated fault segments linked up, the harder it is to identify all the different segments that make up the fault. Thus, different tiers of segmentation along a fault may develop due to the various stages of growth and history of the fault (Jackson and Rotevatn, 2013). As such, the lower tier of segmentation associated with recent surface-rupturing fault offsets (observable along surface fault traces) may not necessarily mimic the higher tier of segmentation associated with long-term fault growth and is observable in the alongfault distribution of cumulative fault offset.
For this study, our scale of observation encompasses the entire length of the fault and maximum vertical offset of the pre-rift basement, thus providing a robust basis for the tier/ classification of the segments. The geometrical segmentation  (Specht and Rosendahl, 1989;Scholz et al., 2020;Shillington et al., 2020) and the tendency for the faults to still accumulate more strain. The arrow lines show the path a fault would follow at different stages (I-VII) as it grows and evolve from initiation to termination (lengthening phase (horizontal arrows); isolated fault initiation and displacement accrual phase (vertical arrows); transition between the normal fault growth phases (I-VI); and at stage VII, fault becomes inactive. The global D-L plots were modified from Kim and Sanderson, (2005)  shows six segments (Figure 4; Hodge et al., 2018) and scarp segmentation shows five segments (Figure 4), whereas the depth-to-basement and overall cumulative throw segmentation show three segments ( Figure 6). Since the cumulative throw is most representative of the long-term displacement pattern along a fault, we posit that the subsurface component of displacement should not be ignored in the investigation of fault zone structure and segmentation as it can contrast observations at the surface. We argue that fault segmentation that is based on scarp height and surface fault trace geometry are more representative of the patterns of most recent coseismic surface rupturing events along the fault. As such, segments with overlapping rupture patches associated with multiple episodes of coseismic surface rupture may show significant cumulative scarp height which may not necessarily be representative of the longer-term displacement history of the other fault segments.
However, we note that although local footwall height (h') is also part of the subaerial component of fault displacement (Figure 3), the similarity of its long-wavelength segmentation pattern with that of depth-to-basement ( Figure 6) suggest that the h' may preserve some information on the long-term displacement pattern along a fault. In essence, our results suggest that for evolving normal faults such as the BMF, the geometrical complexity observed along their surface fault trace might not always be present at depth.
At the scale of the Makanjira Sub-basin (i.e., western branch of the Malawi Rift), the map of depth-to-basement shows that the hanging wall of the BMF exhibits a much wider and deeper zone of subsidence in the northern and central sections than in the south ( Figure 5B). The zone of transition between both regions is marked by a shallowly buried pre-rift basement. This shallowly buried basement also coincides with the footwall of the Bwanje Fault, an antithetic fault within the hanging wall of the BMF (Figure 2A). Overall, we interpret this "topbasement" topography to represent a partitioning of depocenters along the BMF hanging wall in which the depocenters in the northern/central sections were once separated from those in the south by a basement-high. This basement-high has not been significantly buried since the coalescence of the northern/central and southern depocenters. This interpretation is supported (or at least not contradicted) by the relatively shallow basement and distinct hanging wall electrical stratigraphy of ERT Profile 3 (acquired near the transition zone) which contrasts those of the other profiles.

Evolution of Segmentation and Linkage Along the Bilila-Mtakataka Fault
The along-fault distribution of cumulative throw provides insight into how tectonic processes might have contributed to the present-day architecture of normal faults. Our results show the BMF is characterized by a series of relay zones (defined by throw minima) among which the most prominent are those located between the three composite segments ( Figure 6). The lateral continuity of the fault scarp ( Figure 2A) as well as the continuity of offset along the entire fault ( Figure 6) indicate that the relay zones have been breached and that the composite segments are hardlinked. However, among the two major relay zones, the relay zone connecting the Central and Southern composite segments (Kasinje-Chitsulo relay zone in Figure 4) has the least cumulative throw (Figure 6), suggesting that it has not accommodated significant long-term tectonic strain postlinkage of the composite segments.
At the surface, the base of the BMF escarpment (mountain front) is covered by floodplain deposits of the Bwanje River (axial stream) and incised alluvial fan deposits of transverse streams flowing down from the footwall of the BMF (e.g., Figure 2D). This is consistent with observations in the outcrops of the Quaternary Chipalamawamba Beds nearby in Malombe Graben ( Figure 1C; Van Bocxlaer et al., 2012) where axial stream alluvial deposits overlay a deeper unit of lacustrine origin. This deeper lacustrine unit is possibly associated with infrequent flooding of Lake Malawi which deposited thick clayrich sediments into the southern branches of the Malawi Rift (Dulanya et al., 2014;Dulanya, 2017). Most of our resistivity profiles commonly show a general geoelectric stratigraphy in which a shallow low-to moderate-resistivity layer is underlain by a deeper and relatively thicker low-resistivity layer within the synrift sedimentary cover, consistent with geoelectrical stratigraphy onshore of Lake Malawi in the northern Malawi Rift (Kolawole et al., 2018). We interpret that the ubiquitous deep-seated lowresistivity geoelectric layer within the syn-rift stratigraphy represents clay-rich lacustrine sediments deposited during an episode of flooding of Lake Malawi, and that the topmost low-to moderate-resistivity layer represent floodplain deposits of the axial streams and incised alluvial fan deposits.
In contrast, the local hanging wall geoelectrical stratigraphy at the Kasinje-Chitsulo relay zone (ERT Profile 3; Figures 9A,B) shows a deep sedimentary geoelectric layer that hosts 100-300 Ωm discrete bodies within a background of 32-100 Ω-m, interpreted to be colluvial deposits of less-weathered basement material dispersed within a mountain front alluvial fan stratigraphy ( Figure 2D). The absence of this unique electrical stratigraphy in the other ERT profiles (this study and onshore of northern Lake Malawi) suggests that this relay zone is, perhaps, a major long-lived entry point for alluvial and colluvial sediments sourced from the BMF footwall into the Makanjira Trough subbasin.
Based on these observations, and the collocation of this relay zone with the zone of transition between the northern/central and southern BMF depocenters, we infer that the Central-to-Southern composite segment relay zone is likely the most recently breached relay zone along the BMF. It is possible that after the nucleation of the three composite BMF segments, the northern and central composite segments linked first and continued to localize strain which facilitated the early coalescence of their depocenters. It is also likely that progressive offset on the antithetic faults within the BMF hanging wall (e.g., Bwanje Fault) further amplified the subsidence and broadening of the northern and central BMF depocenters ( Figure 5B). In addition, the results of our study suggest that the Southern Composite Segment has accommodated the least cumulative offset and thus, long-term tectonic strain relative to the other composite segments along the BMF. Thus, we interpret that the BMF segments have been accruing vertical offset over the fault life so far, a portion of the long-term faulting activity has been dominated by the breaching of relay zones and linkage of the fault segments leading to the large-scale lengthening of the BMF.

Implications for Strain Localization During Early-Stage Rifting and Associated Earthquake Hazards
Various studies have investigated the magnitudes of cumulative displacement that a border fault may accommodate before it becomes mechanically incapable of accommodating more strain (Scholz and Contreras, 1998;Goldsworthy and Jackson, 2001;Ebinger, 2005;Corti, 2009;Muirhead et al., 2016Muirhead et al., , 2019. During this period, the border faults accrue throws greater than 1 km and induce extension within their hanging walls, which is then accommodated by the localization of intra-rift faults as slip rates of border faults progressively reduce (Muirhead et al., 2016). However, it has been shown that pre-existing structures may localize tectonic strain along the rift axis during the early stages of rifting prior to significant displacement accrual on the border faults (Wedmore et al., 2020;Kolawole et al., 2021a). The Malawi Rift has been proposed to have propagated southwards (e.g., Scholz et al., 2020); which implies that the border faults in the south (such as the BMF) may be relatively less developed than those in the northern sections of the rift, and that a comparison of the long-term displacement on the northern and central Malawi Rift faults with the BMF may provide a template for understanding the future growth of the BMF.
Further north along the Malawi Rift ( Figure 1B), where the intrarift faults account for a significant portion of the overall tectonic strain , the 170 km-long Livingstone, 140 km-long Usisya, and the 162 km-long South Basin border faults record significant maximum displacements (~7.4,~7.3 and 3.85 km respectively; Figure 11) (Specht and Rosendahl, 1989;Scholz et al., 2020;Shillington et al., 2020). It is possible that these faults may have reached their maximum displacements and might no longer be accommodating increasing strain as seismic slips (Cowie et al., 2005). Nevertheless, we note that the Livingstone, Usisya, and South Basin Border Faults have similar fault lengths as the BMF (Figure 11).
Considering the empirical displacement-length ratio relationships for normal faults (Figure 11; Kim and Sanderson, 2005;Pan et al., 2020), the BMF should have a similar magnitude of maximum cumulative displacement as the northern and central Malawi Rift border faults with which it has comparable lengths, but this is not the case (Figure 11). Based on the integration of both the subsurface and subaerial components of fault offset and dips from ERT and microseismicity cluster, we estimate the maximum displacements along the BMF to be~580-837 m (median of 596 m), which are significantly exceeded by those of other Malawi Rift border faults and global normal faults of equivalent length. Thus, we suggest that the BMF is in fact an under-displaced border fault. Furthermore, we posit that this discrepancy in the displacement-length ratio implies that large magnitude seismic activity along the BMF is possible because there is still room for displacement to be accumulated before the border fault potentially becomes inactive as it progressively attains the displacement-length ratio of the larger-offset northern and central Malawi border faults.
Also, considering the mathematical expression for the fault displacement-length relationship (D = CL n ), the direct proportionality of these two attributes of the fault is dependent on the constant C, which could be influenced by factors such as lithology, rheology, strain rate and direction of tectonic stretching, and strain inhibition or transfer to neighboring pre-existing or subsequent structures (Bergen and Shaw, 2010). In the case of the BMF, we speculate that its underdisplaced history could be a result of strain inhibition by and/or distributed transfer of strain to its synthetic Chirobwe-Ntcheu Fault, the southern sections of the South Basin Border Fault, and the Malombe Fault (prominent fault in the Malombe Graben east of the BMF; Kolawole et al., 2021b; Figure 1C). In addition, we speculate that Malombe Graben could probably be an equally or more active branch of the Malawi Rift and therefore may be accommodating most of the recent strain.
Furthermore, we estimate strain rates along the BMF as a function of cumulative strain (cumulative displacement) against time (initiation of rift faulting at 23 Ma; Van Der Beek et al., 1998;Mortimer et al., 2016;Ojo et al., 2021) using the relationship strain rate = displacement/time. The BMF has lower strain rates (0.025 mm/yr) compared to the most prominent and larger offset border faults along the Malawi Rift which include the Livingstone Fault (0.32 mm/yr), Usisya Fault (0.31 mm/yr), and the South Basin Border Fault (0.20 mm/yr). We suggest that the disparity in displacement may also be due to a combination of lower strain rate of the BMF compared to the Livingstone Fault, Usisya Fault and South Basin Border Fault, the difference in the normal fault growth model, and to inherited heterogeneities of the crust (Williams et al., 2021 preprint).
The presence of potentially recent linkage of composite segments along the BMF indicates a history of prolonged coalescence and offset accrual within the composite segments and a delayed fault lengthening by linkage and coalescence of the composite segments. Thus, we consider a model whereby the BMF could have evolved through alternating stages of lengthening and displacement accrual demonstrating the hybrid fault growth model (Figure 11; Pan et al., 2020) and/or by fault segment linkage Cartwright et al., 1995;Peacock et al., 2017).

CONCLUSION
We integrated SRTM-DEM, electrical resistivity tomography, and aeromagnetic methods to investigate the architecture, cumulative throw distribution, and maximum displacement along the Bilila-Mtakataka Fault (BMF). The results show that although the BMF has attained a possible maximum length similar to the other prominent border faults along the Malawi Frontiers in Earth Science | www.frontiersin.org Rift and global normal fault populations, it has significantly lower cumulative displacement. We provide evidence suggesting "recent" linkage and coalescence of the composite segments of the BMF, indicating that much of the fault life has been spent on the coalescence and offset accrual within the composite segments and a delayed fault lengthening by linkage and coalescence of the composite segments. In essence, the BMF is an under-displaced border fault. Our observations highlight the potential for large magnitude seismic activity along the BMF because there is still room for displacement to be accumulated before the border fault potentially becomes inactive, and it progressively attains the displacement-length ratio of the larger-offset northern and central Malawi border faults.
In a wider context, the multidisciplinary research approach used in this study is useful for investigating the long-term and short-term structural and tectonic evolution of active normal faults in early-stage continental rift settings such as the Malawi Rift. The research approach is especially useful in regions where geological and environmental conditions are unfavorable for the preservation of short-term tectonic indicators [e.g., Rhine Graben and Andean Pre-Cordillera of Western Argentina (Megharoui et al., 2000;Fazzito et al., 2009)].

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.