Consumed from within? Social and ecological drivers of internal clearings in Ethiopian Orthodox church forests

In a largely deforested landscape, thousands of small pockets of indigenous forest remain in Ethiopia's northern highlands. These “church forests” are maintained by Ethiopian Orthodox church communities, and they provide valuable reserves of indigenous species diversity and ecosystem services. However, there are important trade‐offs between the different ecosystem services church forests provide: the size of community‐made clearings within church forests used for prayer and other church services have to be balanced with services derived from intact forest cover and dense vegetation. This research uses a spatial, social–ecological model to examine patterns in church forest size and vegetation density across 2743 church forests in the Amhara Region. Results suggest that larger internal clearings correlate with larger forest area, and clearings also scale with population size. We also find that greater elevation and distance to population centers are associated with higher church forest density, while greater community wealth is associated with reduced forest density and area. These trends initially suggest that internal clearings can support conservation efforts, but they are nuanced by prevalent non‐native Eucalyptus planting and increasing internal deforestation for graveyard space for wealthier communities. Additionally, we find that internal clearings are an important source of omitted variable bias in estimating the relationship between density and area in church forests. As a consequence, incorporating church clearings into social–ecological models leads to meaningful differences in church forest density and area estimates. Ultimately, this study suggests that church forests in Ethiopia may be sustained rather than consumed from within.


| INTRODUCTION
Sacred natural sites, natural areas with spiritual value for peoples or communities, are important for both conservation and the cultural and spiritual values they provide (Berkes, 2009;Ray et al., 2014).
Across the world, sacred sites offer key shelter for plant and animal species (Bossart & Antwi, 2016;Dudley et al., 2009;Mgumia & Oba, 2003), and as a result, the biodiversity of these sites is internationally recognized (Wild et al., 2008).There is evidence to suggest that sacred natural sites in Sub-Saharan Africa may be more effective at habitat conservation than state programs (Strauch et al., 2016), leading some ecologists to advocate prioritizing them for preservation and conservation (Shen & Tan, 2012).Sacred sites also provide benefits to surrounding communities that maintain them, such as water filtration, soil erosion mitigation, and other ecosystem services and economic benefits (Bodin et al., 2006;Hernández-Morcillo et al., 2013;LoTemplio et al., 2017;Millennium Ecosystem Assessment, 2005;Rutte, 2011).Across the largely degraded landscape of Ethiopia's northern highlands, church forests are one such example of sacred natural sites (Wassie, 2002).
Existing literature has examined church forest composition and uses through both ecological and social-ecological lenses.Past ecological studies on church forests have found that environmental and climatic variables are important determinants of vegetation structure and tree species diversity (Aerts et al., 2016;Wassie et al., 2010).Wassie et al. (2010) suggest that elevation and grazing intensity negatively impact forest diversity, while the size of forests has no association with the diversity of species present.Similarly, Aerts et al. (2016) suggest that church forest vegetation composition varies with elevation and rainfall, but there is great and relatively unexplained variation in species composition across forests.In contrast, Berhane et al. (2013) suggest that micro-environmental factors and spatial characteristics are the most important determinants of species composition.
However, many studies acknowledge that church forests are socioecological systems managed by longstanding community-level institutions and suggest that incorporating socioeconomic factors can more clearly illustrate drivers of forest productivity and tree species diversity (Berhane et al., 2013;Reynolds, Collins, et al., 2017;Wassie et al., 2010).Past research suggests that greater livestock presence is negatively associated with both the size and productivity of forests (Reynolds, Collins, et al., 2017;Wassie et al., 2009Wassie et al., , 2010)).Similarly, there is evidence that walls built around church forests result in higher forest regeneration rates, likely the consequence of excluding livestock (Woods et al., 2017).Additionally, Wassie et al. (2010) suggest that wood harvesting is negatively impacting the vegetation structure of wood stands in church forests.However, the existing socioecological literature largely relies upon very small samples of forest patches.
In a nature's contributions to people lens (NCP), church forests are also rare examples of social-ecological systems where the nonreligious, largely material contributions of forests are often secondary or unintentional relative to religious, largely nonmaterial contributions (Díaz et al., 2018).Recent ethnographic studies have argued that church forests are not managed strategically as natural resource systems, but instead are managed to facilitate religious benefits, many of which center on clearings in the forest (Klepeis et al., 2016;Orlowska & Klepeis, 2018).There are two main types of church forest clearings: the central qitsir, which houses the main church buildings and provides congregation service space, and the mehabir (local church associations) clearings outside or bordering the qitsir, which primarily serve as meeting places and graveyards (Kent & Orlowska, 2018;Orlowska & Klepeis, 2018).Additionally, many church forests serve as living areas for monks and nuns (Orlowska & Klepeis, 2018).Religious significance is greatest in the center church building and the surrounding clearings (Kent & Orlowska, 2018); such benefits extend, to a lesser degree, to the broader forest area, the keqitsir wichi or atsed, although the holiness of the forest space is derived from its relation to the church rather than in the sacredness of individual trees (Kent & Orlowska, 2018;Klepeis et al., 2016).A church-centric management scheme hence seemingly favors internal clearings for religious significance over potential ecological benefits arising from greater church forest area, vegetation density, and species diversity (Klepeis et al., 2016;Orlowska & Klepeis, 2018).
Despite the centrality of religious and church benefits, church forests can also provide a myriad of other contributions to church communities.Material contributions can include firewood, charcoal, construction timber, fodder, tools, food, medicine, and tree seedling provision (Amare et al., 2016;Reynolds, Collins, et al., 2017;Sahle et al., 2021), in addition to materials for religious sacraments (utensils, oils, incense, and fruit) (Bongers et al., 2006;Sahle et al., 2021).Nonreligious nonmaterial contributions can include shade, communal gathering space, and honeybee pollination, while regulating contributions include temperature reduction, water filtration, and rainfall distribution (Amare et al., 2016;LoTemplio et al., 2017;Reynolds, Collins, et al., 2017;Sahle et al., 2021;Tura et al., 2016).The types of resources available will depend on the species of trees present, and thus diversity of resources should generally scale with diversity of the forest patch (Reynolds, Collins, et al., 2017).Resource extraction is typically forbidden for non-church functions, although such rules are often poorly monitored and enforced, and differ across churches (Orlowska & Klepeis, 2018;Reynolds, Stave, et al., 2017;Sahle et al., 2021).
Internal clearings vary in size across church forests, and thus help illustrate the outcomes of communities balancing religious contributions with nonreligious contributions.Clearings may look like examples of degradation and deforestation, but they are a necessary aspect of providing religious nonmaterial contributions and enabling churchbased conservation in this context (Klepeis et al., 2016).Several studies have identified these internal clearings as important aspects of church forests' composition (Cardelús et al., 2013;Klepeis et al., 2016;Liang et al., 2016;Orlowska & Klepeis, 2018), and one previous study reported that human disturbances and clearings were negatively associated with forest size, tree species richness, and tree density in a sample of 44 church forests (Cardelús et al., 2019).Additionally, recent studies have raised alarms that church forests may be degraded from outside and within (Klepeis et al., 2016;Reynolds, Collins, et al., 2017;Woods et al., 2017), in particular that increasing wealth may be driving increases in graveyards, headstones for graveyards, and built structures inside the forest (Orlowska & Klepeis, 2018).However, no existing study has quantified the extent and consequences of church forest clearings at a broader scale.
This study takes place on a landscape level within a socialecological systems context, and it is the largest scale study to incorporate clearings into spatial social-ecological models of church forests.
Our study does not investigate church forest management processes explicitly, such as through Ostrom's SES framework (Ostrom, 2009), although existing literature has developed process-driven SES analyses of church forests (Klepeis et al., 2016;Reynolds, Stave, et al., 2017).Instead, we focus on landscape-level indicators and outcomes over a larger number of forests.This scale is critical to understand church forest outcomes across space and inform more in-depth, smaller-N research efforts.
Our study has two objectives: First, we investigate trade-offs in church forest resources, as well as examine the social, economic, and ecological variables associated with these resources.Second, we explore whether or not internal clearings are omitted data points that are consequential for community forestry research, acknowledging that the additional data collection needed to account for clearings in forest GIS analyses is costly.We address these objectives using a sample of 2743 digitized church forests, spatial landscape data, and household data from the Ethiopian Socioeconomic Survey.We model these data with a series of spatial autoregressive (SAR) models to test the associations between social, economic, and ecological variables and ecological church forest outcomes, including internal clearings.was unclear whether the extent had ended, trees within 10 m of any tree included in the forest tracing were included as part of the extent, while trees further than 10 m were excluded.Furthermore, if a church forest was connected to government tree plantations, other preserved areas of forest or local Eucalyptus plantations, the end of the extent was estimated based on differences in tree age and type (with church forest trees appearing as larger trees) and likely church forest shape, informed by an assessment of trees in the local landscape, which helped determine whether tree patches around a forest were the norm or the exception within a given landscape.Each digitized extent was checked by at least one additional researcher to ensure intercoder reliability.Figure 2 illustrates how spatial characteristics can vary greatly across individual forests.
Additionally, we traced the internal clearings of the churches from satellite imagery that most closely matched the traced church forest extents.We define internal clearing as the contiguous space lacking trees that includes and surrounds a church forest's identified church building.This definition includes the central clearing, the qitsir, as well as any mehabir clearings that are connected to it.If the internal clearing connected with the outside of the forest, the tracing followed the edge of the previously traced forest extent.After tracing the forest area and subtracting the internal clearing, we narrowed our subset of analyses to the 2743 church forests that had 0.5 hectares of surrounding forest or greater.

| Tracing adjustment
We adjusted each church forest tracing polygon to account for the slight offsets in the overlap of Google Earth data and other GIS data, potentially a byproduct of the Google Earth mosaic, using Sentinel 2 data as a guide for redelineation.To do so, we adjusted each polygon to maximize the average normalized difference vegetation index (NDVI) of each forest within a 17.5-m radius of the original polygon (Hansen et al., 2013).This process is effective because church forests often have clearly delineated boundaries, because they most often are surrounded by agricultural or grazing lands (Reynolds, Collins, et al., 2017).The effect of this adjustment can be seen in Figure 3.We include variables for distance to market, firewood collection time, household head education, iddir 2 membership, own bicycle, own electric stove, and own television, summarized in Table 1.Distance to market is the distance to the nearest market from the household, in kilometers.Firewood collection time is a measure of how many hours a household spent gathering firewood in the day prior to the survey response.Household head education is the highest education received by the identified household head, separated into a binary "no education" versus "some education" as most household heads had no education.Iddir membership indicates whether or not the household is part of an iddir organization.Iddirs are collective mutual help insurance associations that provide members with assistant in times of emergency or serious shock (Aredo, 2010;Léonard, 2013).
Own bicycle, own electric stove, and own television signify whether or not the respondent owns one or more of the respective assets.Each asset variable is intended as a partial proxy for wealth, while bicycle and electric stove also encompass increased travel mobility and lesser reliance on fuelwood, respectively.

| Landscape data
For each church forest in the sample, we calculated spatial variables in two ways: one that did not account for internal clearings in estimates of the church forest extent and one that did account for clearings (subtracting the cleared area from each church forest polygon).For each approach, we calculated area, edge-to-area ratio (perimeter divided by area, shown in Equation ( 2 2 shows summary statistics for the landscape variables across the church forest sample dataset. We note briefly that there is a slight mismatch in time period boundaries over time (Klepeis et al., 2016).Additionally, the mismatch in survey year is due to staggered enumeration years of the ESS. 3

| Analytical methods
Statistical tests for all analyses were conducted in R Studio 3.4.3using a SAR error model and specifying the W matrix based on a 100 k adjacency matrix.The SAR model accounts for spatial autocorrelation in the error term, which we detected in most of our models via a Moran's I test (Dormann et al., 2007).We selected this model because we think it is most likely that bias in the error term is from unaccounted spatial variables that influence nearby forests in a similar way, such as soil quality, than from spatial lags in the dependent variable (SAR lag), and we do not expect betas to meaningfully vary across space (Geographic Weighted Regression) (Brunsdon et al., 1996;Kissling & Carl, 2008).While it is possible that forest aspects like flora and fauna species are spatially lagged so that forests in close proximity have similar resources, our dependent variables do not capture this finer forest resource data.Additionally, the SAR model is particularly relevant for our ESS subset, because it accounts for correlated error terms within ESS enumeration areas.Finally, we universally apply the SAR error model in order to ensure the comparability of model betas and simplicity of interpretation, even if another form could have a better technical fit for some of our 10 models.We log forest area, edge:area ratio, population density, and internal clearing size in order to obtain normal distributions and address skewed distributions.

| Landscape variables: Full sample
We  we can be more confident of the validity of our bootstrap SAR results.
(2) As a robustness check, we run the landscape SAR for the 231 forests within the bootstrap SAR, to indicate whether or not a change in betas is due to the subsample.We also report the distribution of p values for the landscape variables across bootstrap iterations.(3) Finally, we also note that even with a SE bias, the bootstrap approach is preferable to alternative approaches.Simply applying the mean of each ESS enumeration area to each forest within would create insufficient variation and impose a strong assumption that each household influences each forest in its area equally.Alternatively, omitting ESS data would lead to a notable omitted variable bias in social variables.tively associated with NDVI, although this association is quite small.A moderate 10% increase in clearing size is associated with a roughly 0.5-unit increase in NDVI (scaled to 1000), while a large 50% increase is associated with a roughly 2.4-unit increase.Model 4 shows that forest area is related to edge:area ratio (À), elevation (+), precipitation (À), and clearing size (+).For every 1% increase in clearing size, forest area is about 0.3% greater.Model 5 shows that clearing size is positively associated with forest area, edge:area ratio, elevation, precipitation, and population.Clearing size is also positively associated with NDVI, although this association is negligibly small.Note: Models 1-3 show how results change when accounting for clearings.Model 1 does not account for clearings in estimating NDVI, forest area, and edge:area ratio.Model 2 does account for clearings in estimating these variables but does not include clearing size as an explanatory variable.Model 3 accounts for clearings and adds clearing size as an explanatory variable.Models 4 and 5 apply the Model 3 approach using forest area and clearing size as dependent variables, respectively.N = 2743.*p < 0.05; **p < 0.01; ***p < 0.005.

| Socioeconomic SAR
(NDVI), and internal clearings.Regardless of size, internal clearings and their tie to strong religious value are a major motivating force for church forest conservation (Sahle et al., 2021), if not a definitional necessity of the institution.Our results suggest that clearing size is additionally important, with forests with larger clearings generally being larger.Larger forests may be able to sustain larger clearings and/or larger clearings may enable the conservation of larger forests.
However, even a slight positive association is notable, as a marginal increase in clearing size would necessitate some decrease in area by cutting down trees.mehabir clearings and lump in other disturbances such as weedy and exotic species, while we only trace the center church building clearing.
Our methods limit our ability to represent all clearing disturbances in the forest, but it does allow us to capture the specific effect of central clearings, which encompass most of the clearing space in the sample, and increase our total sample of traced forests.
We also find that church forest area and vegetation density are generally greater at higher elevations, further from market centers, and in lower income areas, while church clearings are generally larger in these same areas.Remoteness from urban centers is generally associated with higher elevations and greater distance to market; we thus interpret the positive association between elevation and church forest area and density as a remoteness proxy, because elevation should typically be a negative predictor of forest area and NDVI due to greater elevation imposing harsher ecological conditions (Reynolds, Collins, et al., 2017;Wassie et al., 2010).Remoteness could also function as a boundary condition, in which distance from population centers serves as a marginal exclusion effect for potential forest users outside the community.In contrast, wealth proxies such as electric stove and tele- suggests that wealth is reducing density through increased pressure for graveyard sites in mehabir clearings, which are only partially captured in our methodology (Orlowska & Klepeis, 2018).In some forests, burial sites may be increasingly located at the forest's periphery in order to reduce clearing expansion (Sahle et al., 2021).
We also find that population density is positively associated with the size of internal clearing, which we interpret as increased demand for space for religious services.This result is consistent with ethnographic evidence that church forest communities tend to prioritize religious services over other forest services, and by extension larger church communities would desire larger clearings.Population density effects are also positive and slight for forest size when controlling for remoteness proxies.This effect is likely because, in higher population areas, church forest expansion over time has been largely through the planting of non-native trees for commercial or construction reasons, most often Eucalyptus trees (Cardelús et al., 2017;Liang et al., 2016;Orlowska & Klepeis, 2018;Reynolds, Collins, et al., 2017).On the one hand, this expansion may indicate that forests may be resilient to population pressures even as population pressures convert other vegetation and scrubland into cropland and housing (Birhanu et al., 2021;Gebrehiwot et al., 2021), which is consistent with evidence that church forest areas have been resilient over time (Ashenafi & Leader-Williams, 2005;Cardelús et al., 2017;Sahle et al., 2021;Scull et al., 2017).On the other hand, Eucalyptus planting can contribute to native forest degradation and thus exhibit a tradeoff between commercial and ecological value that is masked by forest area (Cardelús et al., 2017;Liang et al., 2016).
Beyond the Ethiopian Orthodox church forest context, our study is an example of how internal clearings can enable conservation efforts in a complex social-ecological system.

| Statistical consequences of internal clearings
Our results further show that internal clearings represent an omitted variable bias in the context of Ethiopian Orthodox church forests (Addicott et al., 2022).Our landscape analyses show that both NDVI and forest size are positively related to internal clearing size and that omitting the clearing size introduces a substantial upward bias.In our case, this bias would lead a researcher using Model 1 to make a type I error and conclude a spurious positive correlation between NDVI and forest size.These results are a bit more muddled when we introduce socioeconomic variables.Our socioeconomic subsample loses the significant relationship of NDVI with forest area and internal clearings due to the subsample, although omitted variable bias in coefficient magnitude remains.This omitted variable bias is relevant to existing and future studies, which study the relationship between church forest area and NDVI with larger samples suitable for regression analysis (Reynolds, Collins, et al., 2017).
But are tracing internal clearings worth the sizable time investment?Tracing internal clearings means likely doubling the number of tracings needed for a dataset, but this study suggests that internal clearings are a crucial variable concerning forest management.
Omitted variable bias may be applicable to other sacred groves or spiritual-based forest conservation systems where human use and management leaves a delineated impact on forest layouts.However, internal clearings may not be necessary when using forest area or NDVI as a predictor for non-forest outcomes where there is little reason to think that internal clearings will be related to the dependent variable (Addicott et al., 2022).Omitting clearings may also have little impact on simple summary statistics, because many forests had only a small difference in measuring forest area or NDVI when accounting for internal clearings.The exception is the edge:area ratio calculation, where values notably increased when incorporating internal clearings.

| Limitations
This article is constrained by the limited timeframe of our crosssectional data, leaving the time-variance in church forests rather unclear.While we elucidate associations among forests, we cannot say that forest expansion will causally lead a community to expand the internal clearing or vice versa.Similarly, clearing expansions are clearly not always offset, as evidenced by forests that anecdotally see clearing expansion but no forest growth over time in satellite imagery.
While church forests have been relatively stable or potentially increasing in size since 1938 (Scull et al., 2017) or potentially even 1868 (Nyssen et al., 2009), we should seek to understand whether forests have changed in recent years.Relatedly, it is not yet clear whether church communities preserved existing forests or whether they cultivated forests around existing churches, although the tradition of church forests may have dated back to 400 A.D. (Wassie, 2002).As a consequence, we should interpret forest ecological dynamics in the context of active human management dating potentially to the founding of the church or establishment of the forest itself.Additionally, the quality of our socioeconomic data is limited to the ESS, which only covers a small sample of identified churches and in which we cannot pinpoint the exact distances of households relative to churches.A further limitation is that this study lacks data on soil quality or geological structure, which impacts forest productivity.At the time of writing, the authors know of no available data on these factors for the study area.Finally, our digitization of internal clearings misses any internal clearing outside the central clearing, such as mehabir clearings, as well as any selective cuttings that do not impact the canopy.However, our NDVI estimates partially capture this effect because thinner tree stands produce lower NDVI values.

| CONCLUSION
This study contributes to a growing body of scholarship suggesting that Ethiopian church forests are not necessarily consumed from within; but rather, they are instead sustained from within-with active uses of the forest (including clearings) key to the functioning of the overall social-ecological system.Clearing size is positively associated with forest size in the large sample of church forests studied here, supporting the idea that religious, largely nonmaterial forest contributions through internal clearings may support nonreligious ecosystem service contributions through forest area-rather than reflecting trade-offs in forest management.
From a methodological perspective, this study illustrates how incorporating internal church clearings into ecological models can meaningfully improve the accuracy of those models.The most notable accuracy improvement is the elimination of a substantial upward omitted variable bias in the relationship between forest density and size.
The expanded social-ecological model also suggests that remoteness in the form of higher elevation and greater distance to market centers may lead to marginally greater resource stocks, while wealth may be increasing forest degradation.Population pressures, meanwhile, may increase internal clearing area-but forest area and density appear otherwise resilient.Future studies could investigate this relationship further by tracing the full extent of mehabir clearings within church forests, utilizing higher resolution spatial data with more frequent coverage, establishing a time series of changes to church forest polygons, and incorporating higher resolution sociodemographic data from administering original household surveys in the church communities.
search identified 3185 church forests among the three zones (Figure1).Next, we digitized the extents of the 3185 church forests using satellite imagery between the years of 2003 and 2018 provided by Google Earth Pro.For each church, we digitized the extent of the contiguous forest that surrounds the denoted church building.Barriers such as roads, walls, grazing areas, agricultural fields, and low shrub served as indicators that a church forest's extent had ended.When it F I G U R E 1 Study church forests within the West Gojam, East Gojam, and South Gondar zones.
Four sample church forests in the South Gondar Administrative Zone: (a) Quar, (b) Debre Sena, (c) Wuahir, and (d) Gelawdios.The Quar church forest features a more permeable edge.Debre Sena has a harder boundary than other forests.Of the four forests, Wuahir has the most ragged edge, and the forest also features a larger internal clearing relative to the forest's size.Gelawdios is the largest of the four forests, being almost 70 hectares.These images were drawn from Google Earth Pro.[Colour figure can be viewed at wileyonlinelibrary.com] 2.1.3| Survey data We sourced demographic and socioeconomic variables from the Ethiopian Socioeconomic Survey (ESS).The 2015-2016 (Wave 3) ESS consists of over 5000 household surveys conducted across all regional states of Ethiopia, of which 348 households across 36 enumeration areas fall within the study area (Central Statistical Agency of Ethiopia, 2018).The survey was administered by the Central Statistical Agency of Ethiopia (CSA) with support from the World Bank.Each survey site consists of a sample of households within a 5-km radius of the reported coordinate point, encompassing the enumeration area. 1 Using QGIS 3.4.15,we tied the survey responses to all church forests within a 5-km radius of the surveyed household, defined as any portion of the forest falling within the 5-km radius.The resulting sample includes 231 churches with an average of 12.2 households each.While these data are limited in that it misses households nearby the church sample but outside the enumeration zone, the ESS is still the household survey in rural Ethiopia with the highest data quality and geographic reach.

F
I G U R E 3 Normalized difference vegetation index readings and tracings adjustments for four sample church forests in the South Gondar Administrative Zone: (a) Quar, (b) Debre Sena, (c) Wuahir, and (d) Gelawdios.The Quar and Wuahir church forests have lower vegetation density than Debre Sena and Gelawdios.Additionally, each forest features an unadjusted tracing (gray) and an adjusted tracing (black).These images were created from Sentinel 2 data.[Colour figure can be viewed at wileyonlinelibrary.com] )), forest productivity (measured in NDVI, shown in Equation (1)), as well as elevation, precipitation, and population density.We calculated the mean NDVI value of a forest at a resolution of 10 m from Sentinel 2 data provided by the European Space Agency through Earth Explorer (2018), dated March 3, 2018.We used the Shuttle Radar Topography Mission (SRTM) Void Filled Global digital elevation model (USGS, 2018) provided through Earth Explorer to estimate mean elevation at a 1 arc-second ($30 m) resolution for each church forest.We used the Famine Early Warning Systems Network's (FEWS NET) Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) dataset (2018) to calculate the mean monthly precipitation at a 4-km resolution for each church for February 2017 through February 2018, prior to the March 2018 NDVI estimation.NDVI, elevation, and precipitation data for each church forest were calculated using the Zonal Statistics 0.1 plugin for QGIS.Additionally, we drew population density data from WorldPop (2020) at a 1-km resolution.Table apply the SAR model to the full sample of forests (n = 2743), predicting NDVI based on forest area, edge:area ratio, elevation, precipitation, and population density.We run three models: in Model (1), internal clearings are not accounted for when measuring NDVI, forest area, and edge:area ratio; in Model (2), internal clearings are accounted for when measuring those variables but internal clearing size is not included as an independent; and in Model (3), internal clearings are accounted for and internal clearing size is included as an additional independent variable.Differences in variable values between Model 1 (without clearings) and Models 2-3 (with clearings) are shown in , and clearing size by landscape variables across the study area.Models 1-3 vary the incorporation of clearing data in estimating vegetation density (NDVI), with Akaike Information Criterion (AIC) decreasing as more clearing data are included in Models 2 and 3, indicating that accounting for clearings improves the fit of the models.The magnitude of the negative effect of edge:area on NDVI increases across Models 1-3, even despite Models 2 and 3 having a roughly 40% mean increase in edge:area values due to clearing incorporation.Across Models 1-3, elevation has a relatively consistent positive association with NDVI regardless of clearings: from Models 1 to 3, accounting for clearings in forest area and NDVI calculations leads to a roughly 25-point increase in the estimated effect size of precipitation on NDVI.From Models 2 to 3, the added clearing size variable leads to a roughly 7.5-point increase in the magnitude of the edge:area effect on NDVI.Forest area (+) and population density (À) are significant in Model 1 but lose significance when accounting for internal clearings in Models 2 and 3.In Model 3, clearing size is significantly posi- Models 6-8 vary the incorporation of clearing data in estimating NDVI.We see a similar change in landscape variables across Models 6-8 as we see in Models 1-3: as internal clearing data are incorporated, the forest area coefficient is pushed downwards, the population density coefficient is pushed upwards, and the edge:area ratio effect increases in magnitude.However, the robustness of any relationships between NDVI and forest size and/or internal clearings falls off in this forest subsample.Population density also loses robustness in predicting NDVI in Models 7 and 8, mirroring the landscape regressions.For the socioeconomic variables, the influences of electric stove ownership (À), television ownership (À), and distance to market (À) on NDVI are consistently negative across Models 6-8.Consistent with the landscape SAR, clearing size is positively associated with forest size, although it is no longer positively associated with NDVI (failing to pass at least two robustness tests).It is possible that the NDVI and clearing association is lost due to the forest subsample, as it is also not significant in the landscape model of the subsample (Appendix).Landscape variables significant in the landscape Models 4-5 remain significant and consistent in magnitude and direction in Models 9-10.Two exceptions are precipitation reversing sign in Model 9 and elevation losing significance in Model 10.In Model 9, distance to market (+) exhibits a small positive association with forest size.
These results stand in contrast toCardelús et al.'s (2019) study of human disturbance of church forests, which find a negative relationship between disturbance and forest density.Aside from sample size, the difference in results could potentially be ascribed to a different quantification of internal clearings.Cardelús et al. (2019) include all Summary statistics for the socioeconomic variables from the 2015-2016 Ethiopian Socioeconomic Survey, Wave 3 (CSA, 2018) that fall within our sample set.Summary statistics for the spatial variables of study church forests.
to synchronize precipitation with the NDVI year.A coarser resolution would mask the presence of internal church clearings.The church forest data range vary mostly due to spotty satellite coverage of Ethiopia within Google Earth (the highest spatial resolution image openly available during the data generation period), especially when considering cloud cover.This inconsistency is attenuated by the relative stability of church forest T A B L E 1

Table 2 .
Comparing these models helps to demonstrate whether there is explanatory power in incorporating church forest internal clearings.
2.2.2 | Social variables: ESS subsetFor the subset of forests (n = 231) that are within 5 km of an ESS enumeration area, we run a SAR model based on a bootstrap sample.The point of this model is to test for the influence of socioeconomic variables within a limited subset of forests.For each forest, we take a quence, we omit forests that overlap with these five enumeration areas and no other areas but keep those five enumeration areas' surveys for forests that overlap with one of them and another enumeration area.We tested including the omitted forests and found them to notably push SEs downwards across a number of social variables.An additional methodological constraint is that the bootstrap sample only resamples ESS variables, because the landscape variables are unique to each forest.The landscape variable betas and SEs will differ across iterations only according to how much they covary with the shuffled ESS variables.As a consequence, the SEs of the correlations between landscape variable may be downwards biased.We address this bias in three ways.(1) Our earlier landscape SAR model provides a baseline for comparison with our bootstrap SAR results.To the extent that our bootstrap SAR results are consistent with the landscape SAR results,

Table 3
displays the result of SAR models predicting NDVI, forest size Table 4 reports the results of the bootstrap SAR analysis that esti- mates relationships between ecological and socioeconomic variables and church forest NDVI, forest area, and internal clearing size.In Appendix, we assess potential bias in the SEs of landscape variables.Robust significant results are left as is, estimates that passed at least two assessments are marked with †, and estimates that passed less than two are marked with ‡.
Spatial autoregressive error models predicting forest normalized difference vegetation index (NDVI), forest size, and clearing size as a function of landscape variables.
4 | DISCUSSION 4.1 | Resource trade-offs in church forest management In trade-offs between forest size, density, and clearing size, how do church communities maximize forest contributions?Our results do not show evidence of trade-offs between forest area, density T A B L E 3 vision ownership are all associated with less dense forests.It is worth noting that electric stoves could have the potential to offset firewood pressures on church forest tree stands, but we cannot distinguish this potential effect from the apparent wealth effect.Existing work further T A B L E 4 Spatial autoregressive error models predicting forest normalized difference vegetation index (NDVI), forest size, and clearing size as a function of landscape variables and socioeconomic variables from the Ethiopian Socioeconomic Survey.Note: Model 6 does not account for clearings in estimating NDVI, forest area, and edge:area ratio.Model 7 does account for clearings in estimating these variables but does not include clearing size as an explanatory variable.Model 8 accounts for clearings and adds clearing size as an explanatory variable.Models 9 and 10 apply the Model 8 approach using forest area and clearing size as dependent variables, respectively.NDVI is scaled by Â10 instead of Â1000 in Models 9 and 10 to display the smaller correlations.N = 231 per iteration, 200 repetitions.p Values = two-tailed test of coefficient pooled across repetitions.SEs of the beta means across repetitions are shown in parentheses.