How do large wildfires impact sediment redistribution over multiple decades?

Wildfires have become an increasing threat for Mediterranean ecosystems, due to increasing climate change‐induced wildfire activity and changing land management practices. In addition to the initial risk, wildfires can alter the soil in various ways—depending on fire severity—and cause enhanced post‐fire erosion. Usually, post‐fire erosion studies focus on a short time window and lack the attention for sediment dynamics at larger spatial scales. Yet, these large spatial and temporal scales are fundamental for a better understanding of long‐term destructive effects of multiple recurring wildfires on post‐fire erosion processes and catchment sediment dynamics. In this study the landscape evolution model LAPSUS was used to simulate erosion and deposition in the 404 km2 Águeda catchment in north‐central Portugal over a 41‐year (1979–2020) timespan, including eight wildfires each burning >1000 ha. To include variation in fire severity and its impact on the soil, four burn severity classes, represented by the difference normalized burn ratio (dNBR), were parameterized. Although model calibration was difficult due to lack of spatial and temporal measured data, the results show that long‐term post‐fire net erosion rates were significantly higher in the wildfire scenarios (5.95 ton ha−1 yr−1) compared to those of a non‐wildfire scenario (0.58 ton ha−1 yr−1). Furthermore, erosion values increased with burn severity and multiple wildfires increased the overall catchment sediment build‐up. Simulated erosion patterns showed great spatial variability, with large deposition and erosion rates inside streams. This variability made it difficult to identify land uses that were most sensitive for post‐fire erosion, because some land uses were located in more erosion‐sensitive areas (e.g. streams, gullies) or were more affected by high burn severity levels than others. Despite these limitations, LAPSUS performed well on addressing spatial sediment processes and can contribute to pre‐fire management strategies, by identifying locations at risk of post‐fire erosion.

Wildfires are a common threat in Portugal. Mateus and Fernandes (2014) estimated that the accumulated burned area was 4.2 Â 10 6 ha between 1975 and 2012 ($45% of Portugal's land surface), the highest of all southern European countries. Also, Portuguese wildfires had the second largest mean fire size compared to other southern European countries (24.5 ha between 2000 and 2011). Besides an increasing frequency of weather conditions conducive for wildfires (Calheiros et al., 2020), anthropogenic activity such as land use change, induced by socio-economic change and urbanization, also makes an important contribution to wildfire occurrence (Pausas et al., 2008;Shakesby, 2011). In the 20th century in Portugal, part of the agricultural land and shrublands was converted first into Maritime Pine plantations, and towards the end of the century into Eucalypt plantations (Hawtree et al., 2015;Jones et al., 2011;Vasconcelos Ferreira et al., 2010). Initially, the objective of this conversion was to increase the provision of hydrological services such as erosion protection and flood mitigation, however the runoff and erosion increase caused by wildfire disturbances can, in the long run, negate the intended benefits (Carvalho-Santos et al., 2019;Nunes, Naranjo Quintanilla, et al., 2018).
Several studies in north-central Portugal have shown an increase in erosion in burnt areas (Campo et al., 2006;Ferreira et al., 2008;Hosseini et al., 2016;Hyde et al., 2007;Shakesby et al., 1993Shakesby et al., , 1996. However, these studies mainly focus on short-term processes during a few post-fire years, and on smaller plot or hillslope scales (Shakesby, 2011;Shakesby & Doerr, 2006). Nevertheless, several studies already formulate the need for investigation of the long-term effects, including a consecutive number of wildfires (e.g. McGuire & Youberg, 2019), which could provide insights into the relationship between soil degradation, vegetation change, and other landscape processes that generally occur over a much wider time window. In particular, a better understanding is needed of the importance of occasional severe post-fire erosion events in comparison with longterm background erosion processes (Nunes, Doerr, et al., 2018;Shakesby, 2011;Shakesby & Doerr, 2006). Furthermore, investigating a larger catchment scale can help identify locations with post-fire erosion risk to help delineate intervention strategies, as well as examine the transport pathways of sediments from burnt areas which can have important negative impacts on water quality (Nunes, Bernard-Jannin, et al., 2018Shakesby, 2011).
Due to the difficulty of conducting field assessments, numerical modelling studies focusing on the Iberian Peninsula have been conducted to investigate post-fire erosion response. However, these have focused on rather small spatial (patch, hillslope, or small catchment) and temporal (the first few post-fire years) scales (Fernández et al., 2010;Hosseini et al., 2018;Soto & Díaz-Fierros, 1998;Vieira et al., 2015Vieira et al., , 2018Zema et al., 2020). Models used include the empirical Revised Universal Soil Loss Equation (RUSLE; Fernández & Vega, 2016;Fernández et al., 2010), the semi-empirical Morgan-Morgan-Finney (MMF; Fernández et al., 2010) model, and the physically based Pan-European Soil Erosion Risk Assessment (PESERA; Esteves et al., 2012), Water Erosion Prediction Project (WEPP; Soto & Díaz-Fierros, 1998), and Limburg Soil Erosion Model (OpenLISEM; . Esteves et al. (2012) applied the PESERA model to two headwater catchments (9.7 and 100 ha) for a 50-year simulation period; however, the model only simulated on-site erosion, not accounting for off-site effects caused by sediment transport and deposition. Further studies conducted at the headwater catchment scale (1 km) using LandSoil (Pastor et al., 2019), the Soil and Water Assessment Tool (SWAT; Nunes, Bernard-Jannin, et al., 2018, Nunes, Naranjo Quintanilla, et al., 2018, and OpenLISEM  show the importance of wildfires for erosion and sediment yield even at longer time scales (20 and 10 years, respectively). However, the complexity of these modelsand the large number of parameters required to appropriately simulate post-fire impacts-has limited their applicability with sufficient spatial resolution for larger areas, preventing an analysis of the longterm impacts of recurring wildfires for larger landscapes (>100 km 2 ), which are typically the units at which forest planning is made.
Reduced-complexity landscape evolution models (LEMs) are capable of investigating long-term landscape sediment dynamics (Baartman, van Gorp, et al., 2012;Tucker & Hancock, 2010). Thus, LEMs can be applied to investigate long-term (historical) landscape evolution and sediment behaviour in a wildfire-affected catchment under limited data availability conditions. The aim of this study was to investigate how multiple wildfires have affected spatial and temporal erosion and deposition dynamics in theÁgueda catchment (north-central Portugal, 404 km 2 ) using a long-term modelling approach. We applied LEM Landscape Process Modelling at Multi-dimensions and Scales (LAPSUS;Schoorl et al., 2000Schoorl et al., , 2002 to the study area for a 41-year (1979-2020) time period using a wildfire and no-wildfire scenario. In the wildfire scenario, eight major wildfires which occurred in this period were parameterized, with varying severity and spatial extent; while in the no-wildfire scenario, no wildfires were assumed to have taken place.
Results were evaluated in terms of (spatially explicit) erosion and deposition rates over time and how these were affected by (1) land use and (2) multiple wildfire occurrence. autumn and winter months (Hawtree et al., 2015;Nunes, Naranjo Quintanilla, et al., 2018). Mean annual rainfall of 1787 mm was estimated for the period between 1936 and 2010 (Hawtree et al., 2015).
The geology of the Caramulo mountain range is characterized by schist on lower elevations and granites on higher elevated areas (Hawtree et al., 2015;Tavares Wahren et al., 2016). Soil formation in the region generally resulted in shallow soils with low organic matter content and coarse texture properties (Nunes, Naranjo Quintanilla, et al., 2018). Soils are mostly classified as Leptosols between 40 and 75 cm deep, and Cambisols with a depth between 40 and 100 cm (Tavares Wahren et al., 2016;WRB, 2015). Figure 2 shows the land use in 1995 and 2018, with the percentages for each land use type. Land use evolution (where Pine has been changed for Eucalypt between 1995 and 2018) can be seen, indicating the afforestation practices in the late 20th and beginning 21st century (Ferreira, 1997;Hawtree et al., 2015). Recent recorded major wildfires occurred in 1985occurred in , 1986occurred in , 1991occurred in , 1995occurred in , 2005occurred in , 2013occurred in (Hawtree et al., 2015, 2016 and 2017. A large wildfire, therefore, took place every $5 years in the catchment and had a size ranging between 1500 and 9200 ha. Most wildfires occurred during summer (July, August and September), except for the wildfire in 1986 that happened in April.

| LAPSUS model description
LEM LAPSUS is a physically based model that is able to investigate long-term and large-scale spatial landscape evolution (Baartman, van Gorp, et al., 2012;Schoorl et al., 2000Schoorl et al., , 2002. The original model includes topography, soil depth, changing land use, climate variability (annual rainfall), and soil and vegetation characteristics (annual quantities as well as spatial evapotranspiration and infiltration) (e.g. ; in this study spatial burn severity was also included. Considering the multi-decade timespan of this study, weathering of parent material was not included (Alexander, 1985). Figure 3 shows an overview of the model procedure. Model output consists of annual maps of erosion, deposition, soil depth, and runoff, from which temporal and spatial sediment yield was derived.
LAPSUS is a cellular automata model. Routing of runoff and sediment influx towards neighbouring cells is determined using a multiple flow algorithm based on Holmgren (1994): where f i is the fraction of runoff out of a cell going in direction i (equal to the gradient of the slope Λ in the direction i, powered by the p convergence factor) divided by the total sum of slopes of all lower elevated neighbouring cells j, powered by the factor p.
Incoming water flow Q (m 2 yr À1 ) for a particular cell equals the incoming water flow (from upslope cells) plus the rainfall on the cell.
From this, infiltration and evaporation are subtracted and the remainder is used in the calculation for erosion and deposition and transported to the next (downslope) cell.
For sediment transport, the continuity of transport and conservation of mass principles apply. Equal bulk density of eroded and deposited material is assumed. For each transition from cell to cell along length dx (m) of a finite element, sediment transport capacity C (m 2 yr À1 ) is calculated (Equation (2)) as a function of fractional discharge Q and slope tangent Λ (Kirkby, 1971): with discharge exponent m and slope exponent n, and γ a constant for unit conversion with value 1. This is the well-known stream-power equation (Lague, 2014). Transport capacity is compared to the incoming amount of sediment in transport S 0 (m 2 yr À1 ) to calculate the amount of sediment S (m 2 yr À1 ) that will be transported: with dx the cell size (m). Term h (m) refers to the transport capacity divided by the detachment capacity (D, m yr À1 ) in case of erosion, and transport capacity divided by settlement capacity (T, m yr À1 ) in case of sedimentation: where K es (m À1 ) is an aggregated surface factor representing the erodibility of the surface, while P es (m À1 ) is a similar factor for sedimentation potential. The potential for erosion or deposition depends on lumped parameters K es and P es , aggregating different landscape characteristics (vegetation, soil properties, etc.), which are parameterized for different land uses and in this case for burnt areas, both immediately after the wildfire and for the recovery period.

| Data collection and model parameterization
The main input parameters needed by LAPSUS include a digital eleva- 2018; all years given in Supplementary Material S1). These land use maps were adapted to use a less detailed classification and, for the wildfire scenario, to include the effect of the wildfire in the first, second, and third post-fire years based on burn severity. Only major wildfires >1000 ha were included, which covered $80% of the total burnt area in the 41-year period considered. These wildfires happened in 1985, 1986, 1991, 1995, 2005). Burn severity was estimated by deriving the normalized burn ratio (NBR) using Landsat satellite imagery (ESPA-USGS, 2020), which was used to calculate the different normalized burn ratio (dNBR) by subtracting the pre-fire NBR image from the post-fire NBR image: 2.3.1 | Erodibility and sediment potential parameterization The dNBR maps were classified to derive areas of different burn severity (Table 1) and combined with the land use maps. For each land use-burn severity combination, parameterization of lumped K es and P es parameters was based on the USLE vegetation cover factor (USLE_C; Table 2), as vegetation cover is assumed to affect erodibility (i.e. more vegetation cover leads to lower erodibility) and sedimentation potential (i.e. more vegetation cover is assumed to increase sedimentation of transported material). Since LAPSUS does not account for parameterization of the USLE_C value directly, USLE_C values were used as a proxy for fractional change in K es and P es values which are used in the model (Equations 4 and 5). Thus, the absolute value of the USLE_C factor was not used, but the relative differences in USLE_C factors between land use-burn severity combinations ( To derive the total annual time series of infiltration and actual evapotranspiration, the catchment outlet streamflow was used. One of the reasons to focus on theÁgueda catchment is the availability of the (quite unique) historical daily streamflow data , recorded at the streamflow gauging point 'Ponte deÁgueda' at the outlet of the catchment (Hawtree et al., 2015). To derive actual evapotranspiration, the following equation was used: where ET actual is the annual actual evapotranspiration (mm), P is the annual rainfall (m), and Q is the annual streamflow data at the gauging station (mm); deep percolation losses were assumed negligible due to the relatively impermeable geology.
During the dry months (June to September) there were impoundments in the streamflow data due to the closure of a sluice in the river. Missing streamflow or suspicious data from 2004 until 2020 was corrected based on powered regression of streamflow and rainfall data (R 2 = 0.91), based on the whole rainfall and streamflow dataset . By using this regression curve, known rainfall values were translated into predicted streamflow values.
To estimate time series of annual infiltration, the concept of baseflow (mm/yr) was used, defined as the fraction of the corrected streamflow data that is not directly generated after a rainfall event in the form of runoff. As this is the lateral flow through the soil column and not the runoff transported over the surface, it therefore serves as a proxy for infiltration. Estimation of baseflow was done by the recursive digital filter method (Arnold et al., 1995;Nathan & McMahon, 1990). This method is considered robust no matter the stream variabilities or catchment size (Nathan & McMahon, 1990).
Apart from including time series data, spatial variability of infiltration and evapotranspiration was included. Differences in actual evapotranspiration (ET A ) values (mm) (shown in Table 3, values without brackets) between land uses were based on multiplication by the potential evapotranspiration (PET), equal for the whole catchment, and the crop coefficient (K c ), different for each land use. For calculation of the PET, the Hargreaves method was used, using monthly temperature variation plus altitude and latitude (Hargreaves & Samani, 1985). For estimation of K c , an approach from Maselli et al. (2014) was used, that includes seasonality of vegetation cover and soil moisture to account for dry summer months with lower K c values (see Supplementary Material S2 for an explanatory description about spatial ET).
To estimate spatial variable infiltration values between land uses, the inverse curve number (CN) (Boonstra, 1994; USDA, 2019) was estimated (shown as 'Infiltration' in Table 3, values without brackets).
Inverse CNs were determined based on: land use classes (COS classes); soil hydrological groups (4 classes), slope classes (5 classes); the antecedent soil moisture condition (3 classes) and the hydrological condition related to management (3 classes) (Boonstra, 1994). For further explanation of the use of the CN and the parameterization of spatial infiltration, see Supplementary Material S2.

| Model calibration
To calibrate LAPSUS, the lumped K es and P es factors (Equations 4 and 5) for the whole catchment and the m, n, and p parameters (Equations 1 and 2) were adjusted. Unfortunately, the measured  Shakesby et al. (1996) in the north-west of theÁgueda catchment (see Shakesby et al., 1996 for coordinates) and Ferreira (1997) close to Macieira de Alcôba were used to calibrate erosion rates for  were adapted to 1.35, 2, and 6, respectively.   The net erosion simulated by LAPSUS for the burnt scenario was 5.95 ton ha À1 yr À1 on average, whereas for the unburnt scenario the overall net erosion is one order of magnitude lower: 0.58 ton ha À1 yr À1 . This is mostly due to much higher values of total erosion for the burnt scenario: 16.26 ton ha À1 yr À1 for burnt vs 3.46 ton ha À1 yr À1 for unburnt. There was also an increase in deposition:

| Effects of wildfires on erosion and sediment deposition
10.31 ton ha À1 yr À1 for burnt vs 2.88 ton ha À1 yr À1 for unburnt, meaning that a large part of the additional erosion due to wildfire occurrence redeposited downstream, in non-burnt areas. However, the sediment delivery ratio (SDR = net erosion/total erosion) for the burnt scenario was much higher (36.6%) than for the unburnt scenario (16.8%). Thus, more sediment reaches a stream and finally the outlet in a wildfire-affected catchment. Model results suggest that the main erosion driver in this region was wildfire occurrence, with rainfall variability as secondary driver.

| Spatial patterns of erosion and deposition
For the burnt scenario the land use with the largest contribution to erosion in the catchment was 'Eucalypt' (36.9%; 1.78 ton ha À1 yr À1 ),   (Table 5). This variability is particularly high for 'Other Broadleaf', again determined by their presence in the riparian area. Wildfire recurrency also had a large impact on simulated erosion rates. The simulated median erosion rate in almost all land uses was one order of magnitude larger in areas that burnt one or two times compared to areas that did not burn (Table 5). For three or four times burnt areas, the median erosion rates might even be two orders of magnitude higher. For deposition, median values were almost zero, as it was concentrated in small areas. Mean deposition increased with burn frequency for 'Other Broadleaf (Riparian)' and urban areas (i.e. infrastructure), and to a lesser degree for 'Agriculture' until a limit of two times burnt. Note that there was a large variation in the surface area affected by a certain fire recurrence and its related characteristics (e.g. slope steepness). and is a problem in many long-term modelling studies (Batista et al., 2019). It is also common in post-fire erosion modelling (Lopes et al., 2021), where the unpredictable occurrence and short nature of wildfire disturbances limits data collection immediately after the wild- representation make the LAPSUS model results more suitable to show the spatial patterns of erosion and deposition, and give an indication of the order of magnitude of their rates, than to make precise predictions.

| DISCUSSION
Despite these limitations, almost all calibrated erosion rates were of the same order of magnitude as those observed in theÁgueda catchment at the plot (Ferreira, 1997;Shakesby et al., 1994) and small catchment scales (Nunes, Doerr, et al., 2018), as shown in Supplementary Material S3 and Figure 5. Furthermore, our results were within the range of multiple Mediterranean plot studies reviewed by Shakesby (2011), although slightly at the 'high end' of observations: the median first post-fire year erosion rates of $39.2 ton ha À1 yr À1 in this study fall inside the range 1-50 ton ha À1 yr À1 for moderate to high burn-severity wildfires, with a median of 3.7 ton ha À1 yr À1 . Local Moreover, the median long-term erosion rates found in this study for shrublands (5.46 ton ha À1 yr À1 ) correspond with those found in north-west Spain using 137 Cs inventories to assess erosion at burnt shrub-covered hillslopes for a 50-year period, ranging from 6.6 to 6.7 ton ha À1 yr À1 (Menéndez-Duarte et al., 2009). The erosion values found in this study seem therefore to broadly agree with those reported elsewhere in terms of order of magnitude.
It should be noted that the model might also be limited due to the poor representation of several processes in theÁgueda catchment; streambank erosion, soil renewal, and tillage erosion could all be relevant. For streambank erosion, the lack of sediment yield data at the catchment scale prevented the determination of different erosivity parameters for streams in LAPSUS (Baartman, van Gorp, et al., 2012), and existing stone walls in streams were not accounted for. Soil renewal was also not accounted for, but it can be argued that these values are negligible at the 41-year time scale due to the low soil renewal rates in the Mediterranean region (around 0.1 ton ha À1 yr À1 estimated by Alexander, 1985). Finally, tillage in the managed plantation forests of Eucalypt and Pine might be significant due to the short rotation cycles of 10-12 years (Kardell et al., 1986;Shakesby, 2011), estimated at 10 ton ha À1 after each treatment in central Portugal (Govers et al., 1996;Shakesby, 2011); again, in the long-term these values are relatively low compared with water erosion, which concurs with the findings of  for millennial time scales.

| Sediment dynamics at catchment scale
One of the main findings is the high spatial variability of erosion rates, with a high skewness (and a mean one or two orders of magnitude higher than the median) indicating a concentration of erosion in small areas. This can be attributed to different processes operating at different scales, and concurs with the large difference between measurements at the point, plot, or catchment scale mentioned by Shakesby (2011).
The areas of concentrated erosion in theÁgueda catchment often have high runoff accumulation, being therefore subjected to gully formation or, at larger scales, streambank erosion (Figure 6a). This concentration of erosion is common to many large catchments, as described for example by de Vente and Poesen (2005)  (2020), who found a good relation between sediment yield and topographic connectivity modified by burn severity. However, at larger scales, the role of sediment sinks (footslope and floodplain sediment storage) tends to become increasingly important, lowering the sediment yield (de Vente & Poesen, 2005). This has also been proposed for burnt areas by Ferreira et al. (2008), who associated it with a decreasing sediment transport capacity by runoff and streamflow. In There is, however, variability in erosion and sediment deposition within the stream network of theÁgueda catchment (Figure 6), possibly caused by differences in the connectivity of the stream with upslope locations (Borselli et al., 2008;Heckmann & Schwanghart, 2013). Despite this variability, deposition tends to decrease towards the outlet of the catchment, suggesting a sediment build-up in the larger streams. This process could foster 'sediment cascades' that can behave as 'jerky conveyer belts' towards the stream outlet (Cossart et al., 2018;Schoorl et al., 2014), which has also been observed in other burned catchments (Inbar et al., 1998;Keizer et al., 2015;Mayor et al., 2007).

| The significance of post-fire erosion in thé Agueda catchment
The results indicate that wildfires cause long-term erosion and sediment yield in theÁgueda catchment to increase by at least one order of magnitude when compared to natural (i.e. non-disturbed) erosion and sediment yield (Table 5 and Figures 6c and d, 7).
As shown in Figure 8, the relation between rainfall and simulated total erosion was relatively strong, but decreased when wildfires were included in the simulations. This indicates that, when a wildfire occurs, it clearly dominates the effects of rainfall on erosion. This relationship is influenced by wildfire size, and especially burn severity for the fires in 2013 and 2017. When compared with other fires, the larger extent of high burn severity associated with these wildfires corresponds to a larger increase of net erosion, total erosion, and deposition rates for the whole catchment, especially in the first post-fire year (Figure 7).
There were relatively high erosion values even in the second and third post-fire years, which did not occur after smaller and less severe wildfires, which was possibly caused by a slower regeneration of the vegetation or soil (Maia et al., 2012). Interestingly, this difference in sediment processes persists after the wildfire disturbance is negated by vegetation regeneration, as sediment eroded from burnt areas and deposited in the first post-fire years can be available for reentrainment in years without wildfire-enhanced erosion. This process corresponds to earlier research (Inbar et al., 1998;Mayor et al., 2007;Wittenberg & Inbar, 2009), although in theÁgueda headwater catchments Nunes et al. (2020) found that the large rainfall rates can flush these sediments in a few years.
Wildfire-enhanced erosion rates in theÁgueda can lead to soil degradation. Median values in pre-fire years are under a tolerable rate of 1 ton ha À1 yr À1 (Verheijen et al., 2012), but estimated erosion rates for all three years after each fire surpass this threshold by one or two orders of magnitude. Soil depth loss can be especially severe given the generally shallow soils in theÁgueda catchment (Tavares Wahren et al., 2016); roughly 8% of the catchment showed a reduction of 10 cm soil depth after 41 years. Moreover, erosion can lead to a decrease in soil fertility (Bakker et al., 2004) through the selective loss of carbon and nutrients (Hosseini et al., 2017;Serpa et al., 2020). Even though our simulated erosion rates are uncertain due to calibration difficulties (see Section 4.1), this concurs with the generally degraded condition of mountain soils in this study area (Tavares Wahren et al., 2016), many of which are already unsuited for Eucalypt plantation and support less profitable Pine plantations. The effects of repeated wildfires on soil quality have also been observed elsewhere in the Mediterranean (Carreira et al., 1996;Hosseini et al., 2016).
Furthermore, the higher connectivity caused by the wildfires can also have implications for water quality. Ashes and fine sediments produced by wildfires can impact ecosystems and limit human uses, due to associated contaminants such as polycyclic aromatic hydrocarbons, heavy metals, and nutrients (Nunes, Naranjo Quintanilla, et al., 2018;Serpa et al., 2020). While the increase in net erosion can indicate potential problems in the immediate years after wildfires, the increase in background sediment concentrations can indicate persistent contamination problems, especially where related to toxic compounds which can be problematic in small amounts and accumulate over time.

| Contribution of land use to erosion patterns
The results indicate that the most important plantation tree species, Eucalypt and Pine, are the main contributors to the overall sediment budget in theÁgueda catchment (Figure 9a). However, mean erosion rates for natural vegetation (predominantly consisting of shrub vegetation) are much higher than those for other species (Table 5  Areas with other broadleaf species showed much higher spatial variability, with erosion rates being determined by a low number of cells with relatively high erosion values. This can be explained by the location of part of these forests in riparian areas with higher runoff accumulation, and therefore susceptible to channel bank erosion but also to streambed deposition. This spatial variability is therefore a result of the association between these species and channel processes. Finally, 'Agriculture' and 'Urban' had lower erosion rates than forest land uses, even though 'Agriculture' is parameterized as more susceptible to erosion (higher K es in These results complicate a direct assessment of the land use most responsible for erosion based solely on its characteristics, as the relation between both variables at theÁgueda catchment scale seems to be related to wildfire occurrence and recurrence, or their topographic position, rather than local-scale processes such as soil water repellency or soil cover by needle cast (Benito et al., 2003;Shakesby et al., 1994;Walsh et al., 1994). This concurs with the conclusions of Ferreira et al. (2008) that hydrological and sediment connectivity, and their impact by wildfire, drive burnt area erosion response at wider scales.

| Challenges for land management in the future
Increasing global change could further impact ecosystem services and intensify soil degradation and biodiversity loss, which is detrimental to food and water security (Certini, 2005). In fact, post-fire erosion in the region of north-central Portugal is likely to increase, due to intensified rainfall during the winter months associated with a more intense wildfire regime (Carvalho-Santos et al., 2019;Pastor et al., 2019).
Likewise, the wildfires in the last decade (2013, 2016, and 2017) are larger in size and had a larger impact due to higher burn severity. This is partly due to socio-economic drivers such as rural depopulation and (agricultural) land abandonment, leading to encroachment and fostering fuel load build-up, which is perceived to have led to more wildfires (Llovet et al., 2009;Shakesby, 2011).
Areas of natural vegetation were burnt to a larger extent, higher severity, and shorter recurrence times than other forest land uses.
However, burnt natural vegetation areas in the north-east of thé Agueda catchment received less attention in terms of pre-and postfire management. Unlike for natural vegetation, plantations have economic value and are therefore better managed, which fits into the ongoing management paradigm to organize and support productive forests and leave natural areas outside the management structure (Shakesby, 2011). Managing natural areas is, however, important to limit wildfire occurrence, soil erosion, and further degradation. Moreover, better management can also decrease toxic or contaminated ashes that are transported downstream (Nunes, Naranjo Quintanilla, et al., 2018). Hence, managing scrub encroachment in the north-east of theÁgueda catchment can have positive consequences for natural areas, as well as for downstream areas. Approaches can include low investment options such as prescribed burning to lower the fuel load (Khabarov et al., 2016;Shakesby, 2011), or more conventional but expensive options such as grazing or wildfire breaks (Raftoyannis et al., 2014). Similar options can support managing fuel loads in plantations, together with decreasing tree density (Raftoyannis et al., 2014), although plantation densities in theÁgueda are already relatively low. The results of this work can help identify priority management areas and can support post-fire management with emergency post-fire control by mulching Prats et al., 2012), allowing the direction of these measures to areas with greater susceptibility to erosion.

| CONCLUSIONS
To the best of our knowledge, this is the first study to quantify postfire erosion and map sediment transport and deposition patterns at long temporal (41 years) and large spatial (404 km 2 ) scales. Model calibration was considered good, but lack of data-common in postfire erosion assessments-limited a full validation of the LAPSUS model. We found that wildfires significantly increased erosion in thé Agueda catchment in the last 41 years compared to general background erosion processes by at least one order of magnitude. The first post-fire year had a substantial contribution to erosion, and the postfire erosion increased with burn severity and wildfire recurrence.
Wildfires also amplified the background sediment yield in non-postfire years, due to the increased availability of sediment build-up in the catchment. Simulated post-fire erosion rates varied spatially, and showed a large skewness, indicating that erosion is concentrated in well-connected areas, either due to topographic or burn severity patterns. Median erosion values were within the range of those reported in the literature for plot studies, although close to the higher end: 5.95 ton ha À1 yr À1 in the long term, which surpasses a sustainable threshold of 1 ton ha À1 yr À1 , indicating the potential for longterm soil degradation. While there was differentiation of erosion rates between land uses, this was more related to their topographic position and susceptibility to wildfire (both wildfire severity and recurrence) than their intrinsic characteristics, pointing out the importance of sediment connectivity-and in this case, the changes to connectivity caused by wildfire disturbances-for erosion patterns at larger scales. Finally, we think that due to their spatial outputs, reducedcomplexity models such as LAPSUS can function as a tool to identify locations at erosion risk that can be adopted by land managers in pre-fire and post-fire management strategies.

DATA AVAILABILITY STATEMENT
The input and calibration data used in this study is not our own and the source is referred to in the appropriate sections. The output data is presented in maps and tables in the text. We do not intend to share these results further. Interested readers may contact the corresponding author for enquiries.