Examining the role of solar activity, climate, and the socio‐historical context in high all‐cause mortality (northern Portugal, 1700–1880)

The far‐reaching impact of the Sun on Earth's climate and on people's health and well‐being is a poorly understood and non‐consensual scientific issue, with empirical literature stressing the need to expand the knowledge of such relationships. Here, the interplay between solar activity (SA) and climate, and its likely cascading effects on all‐cause mortality, were examined at several time scales. To this end, the parish records of Braga (1700–1880) and Torre de Moncorvo (1700–1850), in two different geographical locations of northern Portugal (Iberia, SW Europe), were used. Crude mortality rate (CMR) and winter–summer ratio (W/S) values were computed to characterize mortality patterns/trends and couple them with potential relevant drivers: total solar irradiance (TSI) as a proxy of SA, the North Atlantic Oscillation (NAO), and key historical events. What emerged, albeit incomplete, was a complex picture of death deeply embedded in people’s physical and socioeconomic environments, at a time when ubiquitous poverty (and co‐morbid malnutrition) was the most inveterate cause of ill health. After identifying the positive mortality episodes in both municipalities, their incidence was found to be higher in periods of weakened SA (normal/grand minima). Standard inference statistics were used to estimate the significance of the observations. The highest CMR peaks matched not only with wars but also with known wide‐ranging mortality crises, which seem to have been triggered by major agricultural production shortfalls, followed by substantial increases in food prices, driven, in turn, by climate deterioration, including extreme weather occurrences. The outcome was social unrest, famines, and outbreaks of infectious diseases, heightening the death toll. The influence of prominent solar/climate variations was investigated using wavelet transform coherence analysis (WTC). The results showed (multi)decadal oscillations in both (TSI and NAO) somehow regulating mortality. But the WTC analysis also estimated SA signals in low‐frequency mortality dynamics disguised by time‐varying determinants, where distinct players of space weather might have been implicated.

The influence of space weather (Baker 1998) and solar activity (SA) on climate variability and change, as well as on coupled natural-human systemsaffecting social, economic, and political structures or biological processes, with implications for mortality (and morbidity)has been the source of a long-standing debate. Over the past two millennia, solar forcing, reinforced by major volcanic eruptions, has shaped several commonly recognized climatic periods; among them, the 500-year Little Ice Age (LIA; c. 1300-1850) that considerably cooled the Earth's climate (Brooke 2018). Such cooling clearly varied over time and space, but it represents a dominant climate event, with important environmental and historical consequences (White 2014). This 'neoglacial' was aligned with, and seems to have been caused by, relatively long solar super-minima (grand minima)the Spörer Minimum (1421-1550), the extended Maunder Minimum (MM;1618-1723, and the Dalton Minimum (DM; 1790-1830)which are part of the c. 2200-year Hallstatt cycle and cluster at its lows (Usoskin 2017;Brooke 2018). These grand minima formed longerperiod envelopes of SA than the well-known quasiperiodic 11-year sunspot (Schwabe) cycle (SC) and roughly coincided with three particularly cold periods across much of the globe (DeGroot 2018). During such cold(est) (sub)periods of the LIA, falling temperatures both caused and reflected multiple changes in atmospheric circulation, namely in the North Atlantic Oscillation (NAO). This, in turn, altered regional and seasonal patterns of precipitation, increasing the likelihood of extreme hydroclimatic events, such as catastrophic droughts or torrential rains (e.g. Brázdil et al. 2010;DeGroot 2018). Researchers have depicted shifts in the mean or in the variability of the regional climate at periods of grand minimanamely in western (and central) Europe, where moist snowy winters, cold springs, cool and rainy (mid-)summers, and less warm anticyclones in autumn prevailed (Pfister et al. 2018)as deeply affecting past societies: either by triggering (or worsening) regional harvest failures and deaths of livestock (due to shortage of pasture), which have conducted to relatively prolonged food supply scarcities, malnutrition, and famines, harming human immunity and increasing vulnerability to disease; or by leading to a concatenation of eventsviolence (war and rebellions), with deaths in action, the accompanying ravaging of civilians, the disruption of farming, the invading armies confiscating or destroying whatever food remained, the outburst of epidemics (typhus and typhoid fever especially, but also dysentery or cholera) in the wake of the devastation. This culminated in mortality crises, when deaths reached quite atypical levels, all over the 'Early Modern' world (Zhang et al. 2007;Camenisch et al. 2016;DeGroot 2018).
Besides starvation and conflictsacting as sources and dissemination of epidemicsmany human diseases, from cardiovascular to respiratory illnesses, along with altered dynamics of infections, are coupled with climate change by impacting pathogens, vectors/hosts or transmission routes (e.g. Patz et al. 2005;Altizer et al. 2006;Wu et al. 2016;Webb 2018;More et al. 2020). Moreover, research undertaken in the field of heliobiology has identified quasi-periodic changes in the Sunenergetic fluxes from the solar and geomagnetic environmentaffecting, both positively and negatively, a wide range of human health and wellness aspects, differentially influencing mortality rates and patterns from various causes (e.g. Davis & Lowell 2006;Babayev & Allahverdiyeva 2007;Mendoza & Peña 2010;Stoupel et al. 2011Stoupel et al. , 2015Cornélissen et al. 2015;Vieira et al. 2018Vieira et al. , 2019Zenchenko & Breus 2021 for a review). Though the heliogeobiological results obtained up to now indicate the nervous and cardiovascular systems as the most clearly affected, the (plausible) underlying biophysical/ biochemical mechanismsfor instance, involving the alteration of the pineal melatonin functions and the Schumann resonancesare far from being resolved (Breus et al. 2016;Krylov 2017), and the scientific issue is highly controversial.
In this work, the overall goal was to investigate the potential influence of SA on the mortality rates of two municipalities of northern Portugal, Braga and Torre de Moncorvo (from now on Moncorvo), over the period 1700-1880, making use of parish data earlier compiled by two historians: H. David (1992) for Braga and V. Tavares (1997a) for Moncorvo. This 180-year period, covering the final phase of the LIA and the transition to the post-LIA warming, was as well an interesting period of time when health care and immunization in Portugal had little or no effect on mortality, which was more dependent on nutritional status and sanitary conditions (David 1992). We focused on high mortality distributions (that is, above-average values; seasonal and yearly), further exploring the following issues: (i) any significant (common) signal between solar/climate-related parameters and mortality in both municipalities throughout the period of records; (ii) any significant difference in the incidence of increased mortality against the background of long-term variations in SA, namely through the quasi-11-year cycle (e.g. minimum vs. maximum) or during episodes of grand minima compared to epochs of normal (cyclical) behaviour; and (iii) the ability (or lack thereof) of temporal variations in mortality to capture some well-known (multi)decadal quasi-periodicities ('cycles') in SA that may potentially modulate terrestrial phenomena. This seems to be a pertinent topic, motivated by the steady amplitude decline of the former quasi-11-year solar cycles 21-24 (solar cycle 24, SC24, was 'the weakest in 100 years'; Svalgaard 2020) and the unusually long and deep recent solar minimum in December 2019 (starting the SC25; NASA 2020; NOAA 2020), raising concerns that the Sun may be heading for the next grand minimum of activity. In that case, it would likely be a driving force for global mean temperatures and regional climates. For some authors, the terrestrial cooling forced by the new grand minimum will not be able to entirely offset the expected global warming from anthropogenic greenhouse gas emissions (e.g. Feulner & Rahmstorf 2010); for others, it will bring to modern times the unique low SA conditions of the MM, with important implications for different parts of the planet on vegetation growth, agriculture, food supplies and heating needs (e.g. Zharkova 2020).

Regional setting
The municipalities of Braga (~183 km 2 ) and Moncorvo (~532 km 2 ), both located in northern Portugal (Fig. 1), can be tracked back to pre-Roman and Medieval (13th century) times, respectively. In the 18th-19th centuries, the Moncorvo region in the northeast of Trás-os-Montes had no industrial or (important) commercial activities, and Torre de Moncorvo was a small village of 1500 inhabitants (Sousa 2018). In the Minho province, Braga had a population of~40 000-50 000 (1795-1878), distributed across the city (urban parishes 14,15,21,22,23,24 and 34;Fig. 1B) and its surrounding countryside. Of the urban parishes, two (23 and 24) occupied almost entirely the city centre, while the others comprised fields and farms to where the city would grow (David 1992). Until the mid-19th century, Braga was a 'rural city' (in the sense that the boundaries between urban and rural were not yet clear), with the local newspaper concerned about people injured by pigs and cows that roamed the main streets (Sá 1993). It was a territory dominated by great landowners (Archbishops of Braga, Benedictine Congregation, Order of Christ) and their economic interests (Durães 2002). The two municipalities had a dominant agrarian population, with a large number of peasants/tenants and daylabourers living close to subsistence, a widespread condition in Portugal, as in the rest of pre-industrial Europe (Wailes 1997).
Regarding climate, continentality, primarily affected by latitude, distance to the sea, and atmospheric circulation, is one of the main features distinguishing the climatic conditions of the two regions (greater continentality in Moncorvo). This, together with the presence of a mountain range running almost parallel to the coastline, isin the words of the Portuguese geographer A. Girão -'what makes Minho the rainiest province in Portugal and Trás-os-Montes a region with an excessively cold and dry climate' (Girão 1941: p. 176). To illustrate regional climate, some available data from the meteorological stations of Braga andMoncorvo (1951-1980 period;INMG 1990INMG , 1991 are provided in Table 1.

Data sources
As mentioned above, the applied monthly-resolved mortality time series were built based on the demo-  (David 1992). B. Braga in the early 21st century (2012) (Carta Administrativa Oficial de Portugal, CAOP2012.1, Direção-Geral do Território, www.dgterritorio.gov.pt/); urban parishes are underlined. C. Torre de Moncorvo and its relationship with the medieval counties of Santa Cruz da Vilariça, Mós and Urros (crossed line = limits of the three medieval counties; dotted line = limits of the current municipality of Torre de Moncorvo; circles with a cross = parishes in the 13th century; Tavares 1997a). D. Torre de Moncorvo in the early 21st century (2012) (https://populacaodistritodebraganca.jimdofree.com/torre-de-moncorvo-1/as-freguesias/).  1951INMG 1990INMG , 1991 graphic works of David (1992) and Tavares (1997a), addressing, respectively, mortality crises in Braga (1700-1880) and Moncorvo (1700-1850). They utilized, as data sources, the records of individual ecclesiastical parishes (a total of 60 in Braga and 17 in Moncorvo), following the methodologies of Dupâquier (1979). References to the age and the direct/specific cause of death were deficient or absent. Both scholars considered years of 'general' crisis (as opposed to 'local') the ones in which at least 25% of parishes with information were affected by crises (i.e. registered an excessive increase in deaths), whatever their intensity. In addition to mortality data, the following time series were also used: (i) the Portuguese national population (Palma et al. 2017;Maddison Project Database, version 2018, (Luterbacher et al. 2002;monthly data). TSI is a common and useful proxy for estimating the SA evolution over long time scales, which corresponds to the spatially and spectrally integrated radiative output of the Sun incident at the top of the Earth's atmosphere (and normalized to a distance of one astronomical unit, 1 AU, from the Sun; Kopp & Shapiro 2021), whereas the NAO is an atmospheric mode that controls a large amount of the winter European climate variability (Deser et al. 2017). A negative (positive) NAO brings more (less) precipitation and lower (higher) average temperatures to southern Europe (Hernández et al. 2020).

Data processing
The parishes with raw mortality data sets exhibiting gaps were not considered for this study. Therefore, the analyses rely on 36 parishes from Braga and 10 parishes from Moncorvo, whose records were aggregated to give monthly totals for mortality in the two municipalities. The crude mortality/death rate (CMR), defined as the mortality rate from all causes of death for a population, was calculated for the monthly data, taking into account the reconstructed national population data sets of the Maddison Project (Palma et al. 2017;Maddison Project Database, version 2018, Bolt et al. 2018  (1) The population of Portugal was used in Equation 1 as a proxy for a 'standard' population, since the size of the resident population for a given year in Braga/Moncorvo over the specified period of time is unknown. Proxy data are quite common in statistical analysis when direct measurements are not available. In this case, we assume that the dynamics of the Portuguese total population is mimicked by the two locales studied here. We acknowledge the limitations derived from this approach when standardizing our data.
Seasonality is critical in understanding climate, so in order to comprehend how seasons (sensu lato) change and influence annual climatic trends, we grouped CMRs at different time scales, which represent broadly cold/hot periods: November to February (NDJF), i.e. the four cool-season months that are expected to affect more closely human health; March to May (MAM); June to August (JJA), including the three hottest/driest summer months; September to November (SON); and also the hydrological year: October to September (ONDJF-MAMJJAS); October to March (ONDJFM), covering the colder half of the year; and April to September (AMJJAS). These CMR values were after standardized by z-score transformation, with a mean of zero and unit variance. Finally, the number of positive mortality episodes, i.e. the z-scores above zero (corresponding to the mortality values greater than the mean), was counted and assembled according to their occurrence over the quasi-11-year SCdivided in line with the expected greater influence of either weak(er) (minimum 'halfcycle') or strong(er) (maximum 'half-cycle') SA (see Fig. 2)or through grand minima (Maunder and Dalton), both as expressed by the TSI data from Coddington et al. (2016). It seems important to point out that the waveform of Fig. 2 (representing two complete SCs; X and X+1) is only a simplified, smoothed representation of the SC that does not take into account the fact that SA has significant variations on time scales shorter than 11 years, namely those that can produce double-peaked maxima (Hathaway 2015 and references therein). The 'split' points (Pi) shown in Fig. 2, corresponding to half of the fall (P1, P3, P5) and rise times (P2, P4), are intended to discern what we consider here as the minimum 'half-cycles' (P1 to P2; P3 to P4) and the maximum 'half-cycles' (P2 to P3; P4 to P5). The reason for the subdivisions of the 'half-cycles'downwards minimum (down_Min) and minimum upwards (Min_up); upwards maximum (up_Max) and maximum downwards (Max_down)was to obtain segments of approximately similar duration and to look for potential differences between them as regards high mortality incidence. The Min_up and down_Min represent each cycle's initial and final active phases, while the minimum (Min) is when the Sun is least active ('quiet Sun'). It is known that the old and new cycle usually overlap in time, making SA at Min a superposition of two consecutive cycles (Leussu et al. 2016). Together, the Min_up and up_Max stages match the typical ascending/rising phase of a quasi-11-year SC, while the Max_down plus down_Min is the descending/declining phase that follows the maximum (Max) activity period. Seasonal fluctuations in mortality were assessed by using the winter-summer ratio (W/S; Rau 2007). This is an index for seasonality and its changes over time, with values of W/S = 1 indicating no differences between deaths in summer and winter, and values of W/S >1 meaning more deaths in winter than in summer (and vice-versa). Following previous works (e.g. Alcoforado et al. 2015), March was included in the winter period (December to March; DJFM) and September in the summer period (June to September; JJAS).
Frequentist (two-tailed) hypothesis testing for single proportions and for two independent proportionsztest; |z| critical value: 1.96; 95% confidence interval (CI): Wilson score with continuity correction (n ≤40) and Agresti-Coull interval (n >40) (Brown et al. 2001) was conducted to estimate the statistical significance for the observed differences: (i) between minimum (down_Min+Min+Min_up = Total_Min) and maximum (up_Max+Max+Max_down = Total_Max) stages of SA in the high mortality (population of CMR z-scores > 0) incidence of Braga/Moncorvo; (ii) the same as in (i), but between normal minimum and grand minimum (Maunder plus Dalton) times of SA; (iii) the same as in (i) and (ii) between the ascending and descending phases of the SC; and (iv) seasonality of deaths in the two municipalities. Expected data under the null hypothesis (H 0 : no difference) were tested at the 95% confidence level. The Mann-Kendall (MK) (non-parametric) test was employed to determine the presence of monotonic trends in the (whole) time series (H 0 : no trend in the series). This was done taking into consideration the presence of autocorrelation and removing its effect, as proposed by Hamed & Rao (1998). During the procedure, the Sen's slope (Q)identifying the trend and showing as well its magnitude (Sen 1968)was also computed. The standard normal homogeneity (SNHT; Alexandersson 1986) and Buishand range (Buishand 1982) (parametric) tests, and Pettitt's (non-parametric) test (Pettitt 1979) were used to distinguish inhomogeneity over time (H 0 : no shift in the data); the good of using all three is that they should complement each other (Arikan & Kahya 2018). These change point detection tests recognize a stepwise shift in the average (a break) and calculate its significance using Monte Carlo resamplings, being also able to pinpoint the year where a break is likely (location-specific tests; Wijngaard et al. 2003). They can generally detect inhomogeneities in all parts of the time series, though the first test (SNHT) is known for working better in the extremesi.e. it discovers breaks more easily at the beginning and at the end of the serieswhile the Buishand and Pettitt tests are more accurate in detecting changes in the middle of the series (Hawkins Fig. 2. Schematic representation of two complete (consecutive) quasi-11-year solar cycles (smoothed line)indicated herein by the total solar irradiance (TSI)and their partitioning according to the expected dominance of the effect of either low(er) levels (minima) or high(er) levels (maxima) of solar activity (SA) over time. This was taken into consideration in the counting (and categorization) of the positive mortality episodes (standardized CMRs greater than zero) in Braga and Torre de Moncorvo, as detailed in the text. The high mortality during the minimum 'half-cycle' is represented by the down_Min, Min and Min_up (Total_Min; 1977; Wijngaard et al. 2003). Moreover, supposedly, the principle of the Pettitt test (distribution-free rank based) makes it less sensitive to outliers than the other two tests (e.g. Wijngaard et al. 2003;Hänsel et al. 2016;Mallakpour & Villarini 2016). Computations were carried out with the software Addinsoft's XLSTAT v. 2021.1. For all tests, a 5% significance level was used.

Wavelet coherence analysis
Wavelet transform coherence (WTC) is a time-frequency analysis method based on continuous wavelet transforms (CWTs), which is quite effective in the examination of transient and nonstationary signals (e.g. Torrence & Compo 1998). The wavelet coherence can be conceptualized as a localized correlation coefficient in the time and frequency space (values ranging between 0 and 1), with WTC allowing the identification of intermittent crosscorrelations between two time series at multiple scales (Grinstedetal. 2004). Accordingly, WTC spectrashowhow the CWTs covary, with arrows designating the relative phase, so that one can recognize periods of covariance with a common phase relationship, suggestive of causality. We applied WTC to examine temporal variability in the relationship between solar and climate-related forcing variables -TSI (Coddington et al. 2016) and NAO (Luterbacher et al. 2002) and all-cause mortality, in the light of potential causal links. It was performed with the wavelet coherence toolbox for MATLAB provided by A. Grinsted and available at http://www.glaciology.net/ wavelet-coherence. All the time series were standardized (zero mean, unit standard deviation) before utilization. Detrending procedures were attemptedlinear and/or polynomial (squared and cubic terms)but since they did not change the significance of SA or NAO, the results shown here include the trends identified in each pair of variables.Throughouttheanalysis,thefollowingwereused: (i) a complex Morlet wavelet (with the default dimensionless frequency/wavenumber ω0 = 6) as the mother wavelet function to the WTC estimation; and (ii) a 95% confidence interval (p-values < 0.05) for statistical testing of the estimated wavelet coherence, according to Grinsted et al. (2004). The parameter ω0 governs the relative time and frequency resolution, with ω0 = 6 providing a good balance between time and frequency localization (Farge 1992;Torrence & Compo 1998;Grinsted et al. 2004). The significance level of the WTC was estimated using a Monte Carlo simulation, generating a total of 1000 surrogate data pairs, with bootstrapping and a first-order autoregressive (AR1) red noise model (Grinsted et al. 2004).

Mortality (CMR and W/S) time series
The mortality patterns (CMRs) for winter (NDJF), spring (MAM), (mid-)spring-summer (AMJJAS), and W/S of Braga and Moncorvo (bringing together the data of 36 and 10 parishes, respectively) are shown in Fig. 3, along with solar irradiance, for the entire period (CMRs and TSI transformed into z-scores). Figure 3A-F displays notable fluctuations in the death rates of both municipalities; the range of z-scores varying between 5.92 (Min-Max: −1.63-4.29; JJA, not shown) and 12.04 (Min-Max: −1.08-10.96; MAM) in Braga, and between 5.20 (Min-Max: −2.44-2.76; ONDJFMAMJJAS, not shown) and 6.52 (Min-Max: −2.64-3.88; AMJJAS) in Moncorvo. A set of 25 most important positive mortality episodes (peaks), 15 in Braga and 10 in Moncorvo, was identified in the (seasonal and yearly) standardized mortality distributions, using the interquartile range (third quartile minus first quartile, Q3-Q1 or IQR; Tukey 1977). They correspond to both first-order (z-score values >Q3+1.5×IQR) and second-order (z-scores >Q3+3×IQR) outliers. The dates and the main characteristics of those peaks (statistical outliers) are provided in Table 2, the latter according to the literature (David 1992;Tavares 1997a, b;Barbosa 2001). Considering any z-score above 3.29 as an outlier (Tabachnick & Fidell 2013), the number of relevant outliers would be six in Braga (1705Braga ( , 1735Braga ( , 1754Braga ( /55, 1775Braga ( , 1809Braga ( and 1872 and two in Moncorvo (1811 and 1831), further supporting the IQR results. The most significant, 1809 (z-score 10.96 in MAM months; extreme/second-order outlier), occurred in Bragathroughout the deepest phase of the DM (Fig. 3C, E)in connection with the French Invasions of Portugal (1807-1811) during the Peninsular War (1807-1814; Table 2); a conflict responsible as well for a first-order outlier in Moncorvo's mortality (1811; Fig. 3B, D). Another two extremes of mortality, also noticed in Braga (1705 and 1872; Table 2), can be seen in Fig. 3: the first (z-score 4.35 in NDJF months; Fig. 3A) within the adverse context of the Spanish Succession War (1701-1714); the second (z-score 4.18 in AMJJAS months; Fig. 3E) related to one of the most severe outbreaks of smallpox and measles, impacting the city and its outskirts (Table 2).
Regarding long-term trends for the full period, the MK test detected a statistically significant decreasing trend in the NDJF (S = −2268; p < 0.05; Q: −3.15×10 −5 ), ONDJFM (S = −2184; p < 0.05; Q: −4.51×10 −5 ), and MAM (S = −2516; p < 0.05; Q: −2.35×10 −5 ) mortality time series, along with an increasing trend in the JJA (S = 3500; p < 0.01; Q: 4.06×10 −5 ) series of Braga, while no trend was found in the mortality of Moncorvo (for an α level of 0.05). Homogeneity analysis showed statistically significant breaks in several time series of both Braga and Moncorvo, denoting an abrupt shift (upward or downward) in the mean before and after a single detected change point (Table 3). This occurred for either the full period or for some of the different sub-periods that were identified after repeating the analysis, in an iterative process that ended when no further inhomogeneity was validated(α = 0.05). Significant common change points in mortality were highlighted statistically by the three tests: towards an increase (seasonal/yearly) near the years 1727 (JJA series of Braga; 1700-1753 sub-period) and 1734 (AMJJAS series of Braga; 1700-1812 subperiod), marking the transition to a time interval with the presence of several (first-and second-order) outliers (e.g. 1754, 1770, 1775 and 1809) in the series of Braga  Table 3)this change point taking place after the heavy death toll of the French Invasions. But there were also sharp differences between the tests in the identified most probable year of change, such as, for instance, 1736, 1753, and 1861 (upward) in Total solar irradiance (TSI) for the same period is also shown (dashed line); CMR and TSI both standardized by z-score transformation. Grouped months: NDJF -November to February, MAM -March to May, AMJJAS -April to September. Symbols: stars show first-orderoutliers (criterion: above Q3+1.5×IQR); triangles identify extremes (second-order outliers; criterion: above Q3+3×IQR), after Tukey (1977). Vertical dotted lines indicate the most likely change points (by date) detected by the homogeneity tests in each time series (up/downward shifts in the average shown by arrows). Table 2. Characteristics of the anomalously high mortality peaks (i.e. first-and second-order outliers; after Tukey 1977) identified in Braga (1700-1880) and Moncorvo (1700-1850) (sources : David 1992;Tavares 1997a, b;Barbosa 2001). Asterisk: this symbol identifies fully epidemic crises. Bold years correspond to the statistical outliers (IQR method; Tukey 1977) further supported by the criterion 'z-score greater than 3.29' (Tabachnick & Fidell 2013 (1753).
Breakdown of agricultural production (1754). High food prices.
Epidemic outbreaks occurring in the late summer would continue through the following months and into early winter (1754).     the JJA mortality series of Braga (Table 3). These change points were likely related to three (first-order) existing outliers, respectively, 1735/36, 1754/55, and 1862 (plus 1864, 1868, 1872/73, and 1876), driven by food shortages following extended droughts and/or epidemic outbreaks, mostly typhus, typhoid fever, and smallpox (Table 2). Additionally, the SNHT found frequency jumps at the very beginning (or end) of some series (full period or sub-periods). This was the case for 1704 in the ONDJFMAMJJAS mortality time series of Braga (1701-1808 sub-period), and in the MAM (1700-1750) and JJA mortality series of Moncorvo (also producing a shift around this date in the AMJJAS and ONDJFMAMJJAS series); and for 1808 in the 1701-1812 sub-period of the ONDJFMAMJJAS mortality series of Braga (Table 3). The most likely reason behind the first date (1704; upward) was the Spanish Succession War and its background, briefly described in Table 2, while the second (1808; upward) had to do with the deadliest year in Braga (1809). It is obvious that the year 1813 was by far the most concordant change point; it was detected in the time series of the two municipalities including the coldest months (i.e. NDJF and ONDJFM; Fig. 3A, B, Table 3). The W/S values (Fig. 3G, H) indicated that mortality in the time frame 1700-1880 was more often higher in winter than in summer, with W/S values >1 during 61% (95% CI: 53.5-67.5%; p < 0.01) of the time in Braga (1700-1880) and for 65% (95% CI: 57.1-72.2%; p < 0.001) of the years in Moncorvo (1700-1850), averaging as follows (meanAE1SD): W/S Braga = 1.13 AE0.33 (95% CI: 1.08-1.18) and W/S TMoncorvo =1.25 AE0.54 (95% CI: 1.16-1.34). Through the period recorded in both series, the number of deaths in Moncorvo was, on average, 25% (95% CI: 18.4-31.9%) higher in winter than in summer, against 18% (95% CI: 11.8-24.2%) in Braga, but this difference was not statistically significant (p > 0.10). A significant linear decreasing trend was detected in both W/S series (Braga: S = −4540; p < 0.0001; Q: −0.003; Moncorvo: S = −1544; p < 0.05; Q: −0.002). Homogeneity analysis showed a significant abrupt shift at the very beginning of the W/S series of Moncorvo (1702; Table 3, Fig. 3H), with the mean level (AE1SD) decreasing from 3.28AE3.28 to 1.24 AE0.45, while several breaks were detected in the W/S data of Braga ( A statistically significant rising trend in TSI was identified for the period under consideration (S = 4899; p < 0.01; Q: 0.002), with the homogeneity analysis showing also significant shifts in the time series, the following dates being the most likely for the (up-ward/downward) change (Table 3) (Table 3).

High mortality vs. quasi-decadal fluctuations of solar activity
The most relevant data comparing high mortalityrepresented by the population of CMR z-scores greater than zeroand SA (indicated by the TSI records) are shown in Table 4. This reveals the number of positive mortality episodes (monthly CMRs grouped by seasons or yearly data), as distributed over the quasi-11-year cycle of TSI, according to Fig. 2. As one can see, the seasons/years with z-scores > 0 did not surpass 45% (95% CI: 38.0-52.5%) in Braga (ONDJFMAMJJAS) and 50% (95% CI: 42.4-58.3%) in Moncorvo (ONDJFM). Considering this higher mortality, statistically significant differences in incidence between minima (Total_Min) and maxima (Total_Max) levels of TSI were noticed, except for the JJA period in Braga, when the difference was not large, and the null hypothesis H 0 could not be rejected. More than 60% of these identified positive mortality episodes -64% (95% CI: 52.8-74.2%) to 72% (95% CI: 61.5-82.1%) in Braga; 69% (95% CI: 57.6-79.9%) to 76% (95% CI: 66.5-86.5%) in Moncorvooccurred when SA was in a relatively weakened condition, during minima 'half-cycles' and grand (Maunder and Dalton) minima (Total_Min; Table 4). It is worth mentioning that, in the same reference period, David (1992) registered 50 mortality crises in Braga, while Barbosa (2001) documented a total of 63 in Moncorvo, of which 68% (Braga; 95% CI: 55.4-80.6%, p < 0.05) and 59% (Moncorvo; 95% CI: 46.9-70.6%, p = 0.208) happened in times of weak(er) SA (Table 4). Analysing, separately, the incidence in periods of grand minima and in epochs of normal minimum conditions (corrected for the respective time frame), higher mortality occurred more frequently in the first (MM+DM), with differences, in general, significant for both municipalities (Table 4), which draws attention to the possible relevance of the (unusual) periods when SA falls below the normal solar minimum levels. Looking at the proportions of positive mortality episodes in the ascending and descending phases of the quasi-11-year SC, shown in Table 4, there is a greater incidence for the descending phase. Most of the Table 4.  WTC results: finding dynamic patterns of correlation between SA/climate-related indices and mortality Figure 4 presents the coherence spectra that resulted from the WTC analyses between the TSI and seasonal (winter, spring and extended winter) mortality, i.e. TSI → NDJF mortality, TSI → MAM mortality and TSI → ONDJFM mortality, in both Braga (left panels) and Moncorvo (right panels). The existence of significant coherence areas at the 95% confidence level in the WTC spectra suggests a causal connection, i.e. SA is likely to have influenced mortality, at the identified scales (frequency/period bands) and specific time intervals (indicated by starting and ending year), over the period under investigation (1700-1880 for Braga and 1700-1850 for Moncorvo). In the panels, one can observe: (i) a shortlived transient significant coherence (black contoured spots) at quasi-decadal and even shorter scales (≤8-16year); (ii) a more persistent one at longer scales, centred at the 32-64-year band (Fig. 4A, E) and at more than 64 years ( Fig. 4B-D, F). Although apparent in WTC spectra, the latter could not be fully resolved, as a result of the length of the mortality time series. In fact, the main feature was a significant TSI → NDJF/ONDJFM (winter/extended winter) mortality coherence in the 32-64-year range of the wavelet scale, spanning almost the whole period of records in Bragaalbeit the significance was constrained by the existence of edge-effects (Fig. 4A, E)which seems to be centred at a longer period band (>64-year) in springtime (MAM; Fig. 4C). This periodic coherence was comparable to the coherent (strong) signal between the TSI and the mortality in Moncorvo from c. 1750 (winter/extended winter; Fig. 4B, F) and c. 1780 (spring; Fig. 4D) onwards, but mostly in a region where the results may not be reliable due to errors associated with edge-effects. These two bands are close to the so-called Double Hale (c. 44-year; e.g. Javaraiah et al. 2005) and lower Gleissberg (50-80-year; e.g. Ogurtsov et al. 2015;Petrovay 2020) 'cycles' of SA, respectively. Overall, the main pattern of coherence revealed TSI variations preceding by far changes in mortality (arrows pointing about 90°down). Significant coherence at the scale of the dominant quasi-11-year SC was highly suppressed during the surveyed period. Nonetheless, in Moncorvo, a fairly sustained significant TSI → NDJF/ ONDJFM mortality coherence, agreeing with the quasidecadal variation of SA, can be observed in two time intervals, c. 1720-1750 (with TSI varying ahead) and c. 1825-1845 (without any lead-lag effect), both showing an anti-phase relationship (analogous to inverse correlation) between the variables (Fig. 4B, F).
The varying NAO → mortality (CMRs) coherence patterns for Braga and Moncorvo (Fig. 5) revealed a prominent anti-phase relationship (arrows pointing left), led by the NAO. Significant NAO → mortality (January) coherence areas in WTC spectra were observed at the 8-16-year frequency band (quasi-decadal scale) in Braga, between c. 1730 and 1815 (Fig. 5A), and at the 16-32year band (quasi-bidecadal scale) in Moncorvo, between c. 1795 and 1825 (Fig. 5B). Analogous coherence was obtained by offsetting the data one month, i.e. by using the NAO for January and Moncorvo's mortality data for February (Fig. 5D). Significant NAO → Braga/Moncorvo mortality coherence was also found for spring season months (Fig. 5C and E-H). In Braga, such an anti-phase relationship was observed in a persistent way at two period bands: 32-64-year, spanning c. 1750/75-1825 (Fig. 5C, E), and 16-32-year between c. 1770 and 1810 (Fig. 5E), extending to c. 1840 (MAM; Fig. 5G). Lastly, Fig. 5F, concerning Moncorvo, points to the prevalence of significant (quasi-bidecadal) anti-phase (April) NAO → mortality coherence in c. 1765-1825, identical to that observed for the MAM season, from c. 1770 to 1820 (Fig. 5H). The key results to be taken from the wavelet analysis are that the NAO probably played a role in mortality during the winter (January/ February) and spring (March/April/MAM) seasons, the two variables exhibiting an inverse correlation, at two long-term scales: quasi-decadal and quasi-bidecadal (widening to the 32-64-year band in Braga). The first for more than 80 years of the 18th century and the first decade of the 19th century (Braga), and the second for variable (shorter) periods (when excluding the region of the spectrum influenced by edge-effects), particularly after c. 1760/70 (Braga and Moncorvo).

Discussion
Higher mortality incidence in relation to external drivers: exploring possible pathways for causality The mortality (CMR) data sets from Braga (1700-1880) and Moncorvo (1700-1850) offered an excellent opportunity to explore solar/climate-human relationships in SW Europe. Firstly, pronounced interannual variations were found in the CMR time series, with a collection of statistical outliers disclosing the existence of seasons/ years in which an exceptionally high number of deaths occurred in both municipalities over the period of interest. This abnormal mortality arose during solar grand minima, namely in association with the three armed conflicts that started in those epochs (Spanish Succession War -MM; French Invasions of Portugal; and Civil War -DM), or related to regional climate instabilityleading to agricultural crises and/or disease outbreaks either at grand minima or in periods of normal SA conditions (Table 2)and giving rise to statistically significant change points (Table 3), as specified in the Results section. These time series also allowed the detection of a significant relationship between increased mortality (CMR z-scores > 0) and reduced levels of SA (normal/grand solar minimum conditions) ( and the y-axis represents the frequency domain (period; in years). The colour scale denotes the magnitude of the coherence power (analogous to correlation: values from 0 to 1). The lighter shaded region indicates the presence of potential edge-effect artifacts. Black line contours designate areas of significant coherence (p < 0.05) against red noise. The black arrows show the relative phase relationship between TSI and mortality: a rightward-pointing arrow represents in-phase coherence (Δφ = 0°); a leftward-pointing arrow denotes anti-phase coherence (Δφ = 180°); a straight down arrow indicates that the first series leads the second by 90°; an upward arrow indicates that the second series leads the first by 90°. that covered 16 quasi-decadal solar cycles. This was achieved by categorizing the CMR peaks (seasonal and yearly) into minimum and maximum 'half-cycles' of TSI as depicted in Fig. 2. Such a hypothetical relationship can be interpreted in the context of the body of research linking SA to either terrestrial climate or human disorders (psycho and physiological). Accordingly, and gathering independent lines of evidence, we suggest two main types of drivers, which will be discussed in the following paragraphs as pathways A and B (both illustrated in Fig. 6).
Pathway A: on the potential influence of heliospheric factors. -The first (pathway A) allows for the impact on individuals of recognized changes in the solar and geomagnetic environment that can act as biological stressors; directly or indirectly, they seem to affect human psychophysiology and behaviour in different ways, by altering regulatory processes, such as melatonin/serotonin balance, blood pressure, immune system, reproductive, neurological, and cardiac system processes (e.g. Alabdulgader et al. 2018 and references therein). Those alterations are likely to become evident over longer time scales, in the form of severe deterioration in welfare and greater mortality (Zenchenko & Breus 2021). The effects of energetic particles (EPP)covering both galactic cosmic rays (GCR), coming from outside the solar system, and solar energetic particles (SEP)are often harmful and can be dangerous to the human organism (e.g. Vaičiulis et al. 2021). Namely, increases in the intensity of GCR-particle flux have been suggested to increase overall deaths (e.g. Stoupel et al. 2011;Vieira et al. 2018). The flux of GCR at 1 AU is modulated by the SC (similar to solar irradiance), with both (GCR and TSI) inversely correlated (Usoskin 2017;Miyahara et al. 2021). Also, the 22-year Hale magnetic polarity solar cyclethe background for the quasi-11-year SC (e.g. Usoskin 2017)affects the modulation of the GCR. Owens et al. (2015) reported that after reversal of the polar field at the solar maximum, the intensity of GCR increases slowly towards the solar minimum in oddnumbered cycles, but rapidly in even-numbered cycles. In line with the previous statements, one possibility for the higher incidence of positive mortality episodes in times of low (minima 'half-cycles') and extremely low SA (Maunder and Dalton conditions) would have to do with enhanced fluxes of GCR (shielded by the solar and heliospheric magnetic field) entering the Earth's atmospherewith the production of secondary cosmic rays (e.g. Usoskin et al. 2009)and presenting an eventual ionizing (high energy) radiation hazard at ground level. Ionizing radiation is much more difficult to avoid than non-ionizing radiation, such as infrared (IR), visible (VIS) or ultraviolet (UV) light, which can easily be shielded out of an environment (Dunbar 2019). Mounting evidence on the biological effects of space weather suggests also that geomagnetic disturbances are an environmental risk factor to humans (Palmer et al. 2006;Dimitrova 2008;Dimitrova et al. 2009;Vieira et al. 2019). The geomagnetic activity (GMA) shows a depen- Fig. 6. Simplified representation of the two proposed pathways by which low/extremely low levels of SA might have contributed to the higher incidence of positive mortality episodes (CMR z-scores >0) in the two municipalities during the surveyed time period. dence on SC, but one that is more complex than that observed in other SA parameters (Hathaway 2015). For instance, GMA tends to remain high during the declining phase of each cycle, which is ascribed to the effect of recurrent high-speed solar wind streams (HSSW) flowing from coronal holes, while in the maximum years of SA, coronal mass ejections (CME) are the major cause of geomagnetic disturbances (Legrand & Simmon 1985;Du 2011;Obridko et al. 2013;Zerbo et al. 2013;Grandin et al. 2019). Therefore, the detectable differences between the rising and declining phases of the quasi-11-year SC (defined by the TSI; Table 4), though borderline significant (p < 0.10), might be explained by the stronger impact over the declining phase (and at the minimum of the SC) of the recurrent GMA usually associated with those specific solar wind events (HSSW). The interactions of this later 11-year cycle GMA with the Earth's magnetosphere are responsible for geomagnetic storms (GMS) on Earth that may be detrimental to human activity and health outcomes (e.g. Gerontidou et al. 2018). Geomagnetic disturbances can worsen existing diseases and appear to be correlated with significant increases in cardiac arrhythmias, cardiovascular disease, and incidence of myocardial infarction related deaths and morbidity (Dimitrova et al. 2009;McCraty et al. 2017 and references therein;Vaičiulis et al. 2021). Besides, GMS have been associated with significant increases in depression, mental disorders and suicide attempts (Kay 1994(Kay , 2004Gordon & Berk, 2003). However, as underscored by Vaičiulis et al. (2021), the biological effects observedwith GMS of different sources (e.g. CME or HSSW) seem to be also very different, the reason why it is essential to take into account their origin for accurate analysis of potential bioeffects due to GMS.
Pathway B: on the potential influence of SA via the AO/ NAO index. -The second pathway (B) is associated with a larger regional scale climate response (colder temperatures) to falling solar irradiance, possible through a forced shift of the Artic Oscillation/North Atlantic Oscillation (AO/NAO) towards the negative index state, and its amplifying mechanisms (e.g. Shindell et al. 2001). The mechanisms proposed so far to explain the projection of solar irradiance variability on regional climate are multi-step and fall into three categories: (i) variations in TSI (VIS/IR) leading to changes in heat input to the lower atmosphere ('bottom-up' -TSI effect); (ii) variations in solar spectral irradiance (SSI), namely in UV radiation, coupled with changes in ozone concentration, affecting the heat budget of the stratosphere ('top-down' -UVeffect); and (iii) variations in the solar wind and flux of EPP (Gray et al. 2010;Lockwood 2012;Thiéblemont & Matthes 2015). A pathway for how solar forcing influences the regional climate of the northwest of Portugal was suggested by the authors in a former work , following the 'dynamical effect' proposed by Muthers et al. (2016), and consistent with the downward transfer of the solar signal from the stratosphere to the surface ('top-down'). It should be stressed that the solar-climate association is an open debate, namely regarding the NAO pattern, with research achieving different results on solar forcing at the quasi-11-year time scale (Gray et al. 2016;Chiodo et al. 2019;Cook et al. 2019). While Gray et al. (2016) detected a tendency for positive (negative) NAO anomalies to follow solar maxima (minima), with a peak response lagged by about 4 years, the other two studies found no consistent temporal link. Hernández et al. (2020) showed that on the decadal time scale a low number (<42) of sunspots was directly correlated with low NAO values, above which this index appeared to be less influenced by SA. On the other hand, Li et al. (2011) argued that on decadal to interannual time scales GMA driven by solar wind is more likely to affect the NAO than solar irradiance, with stronger and more significant correlations earlier described between the two. Their results indicated a statistically significant GMA (aa index)-winter NAO relationship that was negative for small-to-medium aa and positive for medium-to-large aa, and dominated by the aa data from the declining phase of even-numbered solar cycles. The latter not only suggests that recurrent GMA may be responsible for the observed relationship but also links with what we hypothesize here about the higher proportion of positive mortality episodes in those (declining) phases. Given the marginal significance of our results, more (independent) data and other evidence need to be examined to check for consistency.
In this alternative pathway to a logically possible causal relationship, low SA conditions would have prompted unfavourable changes in the regional climate (throughout the AO/NAO), which in turn would have increased the total number of positive mortality episodes (more meaningful for grand minima; Table 4). This might have been achieved directly by physiological reactions to cold, such as increased risks of blood clotting via higher haemoconcentration (cardiovascular and cerebrovascular diseases) and of infections of the airways (respiratory diseases) (Rau 2007), or by means of intermediate factors, like agricultural (and stock farming) production and transmission of infectious agents and diseases (Fig. 6). However, with regard to the causal connections between the AO/NAO, temperature, relative humidity, and winter mortality/cold-related diseases in European countries, the investigations undertaken to date have provided inconsistent outcomes (Messner et al. 2003;McGregor 2005;Almendra et al. 2017 and references therein). Almendra et al. (2017) highlighted that the lack of a significant correlation between the winter NAO and human health (via its influence on thermal stress), as they found in Portugal today (Lisbon area), cannot be extrapolated to other regions of the country or other countries having different geographical/socioeconomic frameworks. Besides the drop of temperature, as a triggering event for winter mortality, several social factors (e.g. housing conditions, exposure to outdoor cold, and inadequate protective clothing) play as well a key role in mediating the annual coldrelated death toll (McKee 1989;Healy 2003;Rau 2007). Another aspect that may help to understand such inconsistency is that, unlike central and northern Europe, the NAO pattern does not exert control in shaping extreme temperatures in mainland Portugal (Santo et al. 2014).
Pathway B: intermediate factors.
-Considering the intermediate factors in Fig. 6, on the one hand, low temperatures and increased humidity have been suggested to facilitate the survival of microorganisms, like bacteria and viruses, in (bio)aerosols (Xie et al. 2006;Lofgren et al. 2007;Prussin et al. 2018). But researchers also emphasized that understanding infectious disease/ infections transmission is a complex problem, as the determinants are multicausal. It is recognized that infectious disease dynamics provide a wide variety of intriguing and unexplained phenomena, with seasonal forcing capable of generating complex multiyear fluctuations in the infection's prevalence, and even chaos (Altizer et al. 2006). Obviously, these findings add uncertainty to our interpretation. On the other hand, the cooling resulting from a (more) negative NAO phase and, perhaps even more decisive, the shortening of the growing seasonimplying insufficient time for the crucial crops (corn, wheat, rye, barley or oats) to reach maturitytogether with increased rainfall throughout it (also induced by the negative NAO), likely disrupted agricultural/stock farming production and the associated revenue. The data in Table 2 allow us to infer that a strong association between harvest failures and high mortality, through a deadly combination of increased food prices, starvation (likely triggering food riots and social unrest), and disease outbreaks, was present at least until 1831 in Moncorvo and 1864 in Braga. The anomalously high mortality peaks identified here (25 in total) all agreed with years of 'general' crises in the two municipalities (David 1992;Tavares 1997a), some of them well matched with recognized longer-term mortality crises in the country and even abroad (Lebrun 1980;Pérez Moreda 1980;Barbosa 2001;Rodrigues et al. 2008;Pfister et al. 2018). Although David (1992) and Tavares (1997a) labelled those crises as 'mixed crises' (in which clearly identifiable epidemics happened in famine years), they seem to be mostly rooted in climate-induced agricultural decline, causing (agro-)economic distress (in both the short-and long-run), as discussed thoroughly byothers (e.g. Zhang et al. 2011;Pei et al. 2013). In such a context, the poor, who represented the majority, were more strongly affected, particularly if failing harvests lasted for a number of years. Undernourishment not only increases vulnerability to health threats (Saunders et al. 2019) but also the risk of hypothermia after prolonged exposure to cold, wet or windy conditions (Cheshire 2016). The anomalous high mortality in the last decades of the MM -1705 (Braga);1709, 1711-1712, and 1716 (Table 2)occurred in a reported period (1680-1720, peaking in the years 1710-1719) of low temperatures and increased cold spellsmore extreme in winter and springin the country (Silva 2019), concomitant with the coldest phase of the LIA in Iberia (Oliva et al. 2018). This (deteriorated) climatic background was exacerbated by the involvement of Portugal in the Spanish Succession War, generating two statistical outliers in mortality (1705, Braga and 1709, Moncorvo; Table 2), and the early change point (near 1704) detected by the SNHT (Fig. 3D, F, Table 3). During the DM, another armed conflict, the French Invasions during the Peninsular War, and its various consequences, was the most striking regarding mortality (Braga and Moncorvo), as previously stated. Yet, consistent with Correia (1938), the typhus outbreak that spread throughout the war claimed far more victims than Napoleonic troops. All this happened in another documented period (1780-1820, peaking in 1810-1819) of more frequent low-temperature anomalies (again more extreme in winter and spring) in Portugal and highly unstable (hydro)climate conditions (Silva 2019), with excessive rainfall leading to insalubrity, from which the armies (in their camps) were also unable to escape (Oliveira 2008). Leaving wartime behind, 1813a significant change year (Table 3)appears to represent a benchmark in the long-term reduction of mortality, particularly in NDJF and MAM (with repercussions in ONDJFM and AMJJAS; Table 3), not only in Braga (where a significant decreasing trend was highlighted using the MK test) but also in Moncorvo (where no trend could be detected). After this date, the steady increasing trend in the summer deaths of Braga and, more specifically, its significant (upward) shift in the early 1860s (JJA, AMJJAS, ONDJFMAMJJAS series; Table 3, Fig. 3E) may point towards an 'adjustment' of the seasonal patterns of mortality to the climatic conditions of the post-LIA phase. Interestingly, in the agrarian historiography of the Entre-Douro-e-Minho region (to which Braga belongs) the year 1813 is referred to as starting a period of slow and progressive recovery of agricultural production indexes, which would continue (with some setbacks) for the following decades up to the 1850s (Oliveira 1979(Oliveira , 1996Silva 2019).
Final remarks on the suggested solar influence pathways. -At this point, it is fairly clear that over the studied period multiple drivers of mortality could have interacted in complex and poorly understood ways at different time scales, generating nonstationarity and nonlinear dynamics, with likely quasi-11-year cycle-tocycle variations. On the relative importance of the two types of proposed influences/(bio)effects on the detected high mortality (pathways A and B), we can only speculate, since obvious problems hinder a more com-plete examination and none can be solved here. Although evidence has been provided that in modern times environmental energetic phenomena (e.g. solar wind, EPP) can affect humans in different ways, depending on their sensitivity, health status, and capacity for selfregulation (e.g. Alabdulgader et al. 2018), such disturbances appear to be quite small in comparison with the archetypal stressors of the 18th-19th centuries. This was a historical period in which reaching adulthood and surviving a few more years continued to be a huge struggle (life expectancy in Portugal would not exceed 35 years in the 1860s; Veiga 2004), due to the succession and range of debilitating diseasesepidemic overlapping endemic; a status quo that persisted in the country to the early 20th century (Almeida 2014)closely accompanied by the poor socioeconomic conditions of its population, almost totally reliant on the regional (national) agricultural production. With this in mind, we may argue that the empirical facts speak in favour of pathway B as the most plausible causal interpretation, although, ultimately, a combination of the two might be suggested. Unfortunately, the paucity of observations in existing death records poses a serious challenge for more thorough research.

Seasonality in mortality over time
The seasonality of mortality (defined by the W/S) showed a pattern of more frequent deaths in winter compared to summer (W/S values >1) dominating more than three fifths of the analysed period. This occurred despite the significant steadily negative trends detected in both W/S series (using the MK test) and the successive downward shifts identified by homogeneity analysis (Table 3, Fig. 3G, H). In Braga, the signs of reversal in the timing of deaths began around 1813 and became more unambiguous in 1845-1853 (to W/S averages <1). Visual inspection of Fig. 3H highlights an important W/ S decrease after 1837 (not validated by any homogeneity test; α = 0.05) in Moncorvo (mean W/S: 0.86AE0.17; 95% CI: 0.77-0.95), at a period in which upward shifts were identified in the average of both the TSI data (1826/27) and the NAO of January (1838) ( Table 3). These observed reverse patterns in the seasonality of mortality could, therefore, be linked to the onset (and evolution) of the recovery from the LIA (Abram et al. 2016;Brönnimann et al. 2019), accompanied by the sharp drop in the frequency of episodes of extreme cold in the wintertime, between 1820 and 1840, in Portugal (Silva 2019). Oliva et al. (2018) described an 1835-1850 slight trend of warming and a gradual staggered increase in mean annual temperature of~1°C after 1850 in Iberian mountains, while Moreno et al. (2016) reported a period of higher spring-summer (MAMJJA) mean maxima temperature (24.1°C), between 1856 and c. 1875, in the Minho region. However, the higher summer temperatures per se would not be the actual cause of people's deaths, the hot weather conditions just providing the basis for the development of certain types of infectious agents; only together with social factors could the disease have spread among humans and wipe out considerable proportions of the population (Rau 2007). Factors such as the enormous unhealthiness of the streets in urban areas, the proximity of stables and dung heaps to housing spaces in rural areas (the animal manure attracting flies, particularly during hot summer months), the poor conservation of food (that was even worse in towns and cities than in the countryside), the consumption of contaminated (stagnant) water, or even the decline in lactation performance and maternal care due to more intense (and heavy) agricultural work, seem to have facilitated the spread of disease during that season of the year (Morgan 1993(Morgan , 2002Leguay 1999;French & Warren 2007;Rau 2007;Rodrigues et al. 2008;Oris & Fariñas 2016). This is also in line with the documented evidence on the timing (JJA period) and nature of the latest anomalous peaks in the mortality patterns of Braga (e.g. 1864 and 1876), pointing to higher infant/child mortality, mainly due to diarrhoeal diseases and common summertime viruses (see Table 2).
A limitation of our analysis is that it does not take into account age profiles within and between municipalities. This was hampered by the preponderance of death records with no statement of age or, instead, one rather ambiguous, especially for children under age seven (e.g. 'angel', 'young age', 'innocent', 'parvulus', with references to stillbirths very rare; David 1992). This information would help to better untangle the drivers of the perceived seasonality in mortality. The 'age group' entails a differential susceptibility and risk of dying from certain infectious diseases (Altizer et al. 2006;Fisman 2012), and it can be as well higly influential on the effectiveness of the human thermoregulatory system (Parsons 2002). Another issue not examined here is the causa mortis, usually also omitted in death records, which was, in fact, generally unknown. This neither allowed us to ascertain the specific cause-of-death nor to address some particularities of the long-term solar/climate influences on mortality (and morbidity) reported in the literature for more recent times.
Forcing factors to mortality long-term variations: the significance of the TSI and the NAO The results of the WTC analyses (shown in Figs 4, 5) indicated that the investigated causal relationships were neither straightforward nor invariable from 1700 to 1880, with coherent time structures and phase behaviour implying an intricate human response (gauged by mortality) to solar/climate forcings. The TSI-mortality coupling (Braga and Moncorvo) reflected mostly antiphase oscillatory synchronization patterns, but with considerable phase differences. Although they suggest a common and persistent long-term solar modulation of mortality (consistent with the frequencies of the c. 44year Double Hale and the 50-80-year lower Gleissberg cycles), the very lagged wavelet coherence (nearly a quarter cycle out-of-phase) indicated factors other than SA as driving the observed CMRs, which might be as well potential sources of cyclicity in mortality dynamics. Being able to distinguish all the potential modifiers or to estimate to what extent each has contributed to the observed time-delayed response would be highly informative. But, even if the necessary data were available, uncovering the causal structure of such a network would be quite demanding, given the noisy background (i.e. substantial year-to-year variability). Likewise, WTC failed to capture the potential effect of SA variability on all-cause mortality on the decadal time scale, probably due to strong high-frequency noise. Still, we may ask whether the observed anti-phase coherence in Moncorvo at the end of the two grand minima (Fig. 4B, F) does not translate, to a certain degree, the human response to the thermal environment (e.g. experiencing lesser and less prolonged cold stress). This would be due to the amelioration of the harsher winters in this geographical location, especially of air temperature, in agreement with the rising trend of the quasi-11-year SC and two identified (upward) change points in the TSI data (1725 and 1826; Table 3). A similar pattern is recognizable in Braga, though weak and transient (Fig. 4A, E), seeming compatible with the arbitration of an external (and widerange) driver of lower mortality, such as strengthened SA, in those periods.
The persistent significant NAO anti-phase coherence with mortality (i.e. more deaths associated with cooler conditions and more rainfall) in winter (January/February) and spring at the quasi-decadal (in Braga) and quasibidecadal (in both Braga and Moncorvo) scales is consistent with the solar Schwabe/Hale cycles. We can interpret these results as suggesting the occurrence of long-term (decadal/bi-decadal) variations in the NAOdriven by two well-known solar quasi-periodicitiesacting as a link between climate (temperature/rainfall) and mortality variabilities. In the absence of reliable historical regional climate data, it is not possible to discriminate between the effects of temperature and precipitation (these can even be opposite when considering different causes of death). Significant quasi-decadal NAO → (January) mortality coherence in Braga, but not in Moncorvo, until c. 1815, close to a statistically significant downward shift detected in NDJF mortality (1813), may probably be due to main differences in regional climatehumid in Braga; cold and dry in Moncorvoin shaping decadal changes in mortality rates. In wintertime, the NAO appears to best explain interannual variations in precipitation, with the influence of other leading atmospheric patterns, such as the Eastern Atlantic (EA), surpassing that of the NAO in the modulation of temperature in SW Europe (Sáenz et al. 2001;Ulbrich et al. 2012). Alternatively, it could be related to NAO-EA interactions, pointing to the EA as a potential source of nonstationary signatures of the NAO (Mellado-Cano et al. 2019). The possibility of a fingerprint of the 22-year Hale cycle in the (winter) NAOmortality relationship during the MM and DM (Fig. 5A, B, D; information on the MM may be distorted by edgeeffects), disappearing soon after, seems to be compatible with the behaviour (amplification and anomalous GCR spikes) of the Hale cycle during periods of extremely weak SA (Miyahara et al. 2009;Kataoka et al. 2012), and expected influences on large-scale climate modes, e.g. inducing negative AO/NAO phases, cascading to human disease and mortality dynamics.

Conclusions
Several data sources and methodologies were combined to explore the role of the physical and socioeconomic environments in explaining mortality patterns/trends (and peaks) in two pre-industrial Portuguese municipalities (1700-1880). The results of the analyses together with the currently available historical testimonies suggested: (i) (multi)decadal oscillations in solar/climaterelated parameters influencing death rates; and (ii) shorter-term detrimental changes in regional climate harmfully impacting agricultural and human (social and biological) systems, which led to increased mortality. Importantly, high all-cause mortality was found to be preferentially associated with periods of weakened SA (quasi-11-year minima "half-cycles" plus grand minima), implying lower levels of radiant energy emitted at all wavelengths, less solar magnetic shielding, and more GCR reaching the Earth. On the face of this, greater frequency of anomalously high mortality over (grand) solar minimum conditions, in Braga and Moncorvo, might have been the aftermath of the cumulative impact of low TSI/high GCR-particle on individuals and their health condition. Two probable causal pathways for such solar influence have been described here. The results of the WTC analysis also highlighted the existence of convoluted relationships between SA, climate, and human mortality (nonlinearity and nonstationarity), concerning the oscillatory modes (periods/frequencies) detected in the analysed data. Additional time series (populational and environmental) and further explorations are needed to support both the accuracy of our inferences/suggestions and their underlying assumptions, leading to a better understanding of these challenging scientific questions. From a historical perspective, major restrictions stand in the way. The under-registration of infant/child deaths (introducing biases in death estimates or in stratified analyses), misleading diagnoses of diseases (preventing the establishment of a parallel with those from more recent times), the advent of vaccination (eradicating many prevalent infectious diseases), and increases in life expectancy (with the emergence of new illnesses associated with ageing and genetic damage) are just a few. To conclude, as the essence of the main subject explored here continues the same, it should not be disregarded by contemporary analysts/policy-makers. Some of their goals might be to make surveillance/monitoring of drivers available (or improve it); formulate and implement strategies to selectively address the problems posed by the Sunclimate interaction on populationsagrarian and nonagrarian livelihoodsand mitigate its demographic consequences (targeting the most vulnerable); and raising (local/global) awareness of these lesser-known natural hazards, particularly among environmental/public health professionals.