longtable

# Integration of biophysical connectivity in the spatial optimization of
coastal ecosystem services111© 2020. This manuscript version is made available
under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-
nd/4.0/

Andres Ospina-Alvarez222Corresponding author: Andrés Ospina-Alvarez, email:
<EMAIL_ADDRESS>address: Spanish Scientific Research Council, Mediterranean
Institute for Advanced Studies (IMEDEA-CSIC/UIB), C/ Miquel Marques 21, CP
07190 Esporles, Balearic Islands, Spain. Mediterranean Institute for Advanced
Studies IMEDEA (UIB-CSIC), C/ Miquel Marques 21, CP 07190 Esporles, Balearic
Islands, Spain.  Silvia de Juan Mohan  Institute of Marine Sciences ICM
(CSIC), Passeig Marítim de la Barceloneta, 37-49, 08003, Barcelona, Spain.
Katrina J. Davis  Department of Zoology, University of Oxford, Zoology
Research and Administration Building, 11a Mansfield Road, Oxford, OX1 3SZ,
United Kingdom  Australian Research Council Centre for Excellence for
Environmental Decisions, The University of Queensland, St Lucia, Queensland,
4072, Australia.  Catherine González  Instituto de Fomento Pesquero (IFOP),
Almte. M. Blanco Encalada 839, Casilla 8-V, Valparaiso, Chile.  Miriam
Fernández  Núcleo Milenio – Centro de Conservación Marina; Estación Costera de
Investigaciones Marinas, Pontificia Universidad Católica de Chile, Alameda
340, C.P. 6513677, Casilla 193, Correo 22, Santiago, Chile.  Sergio A.
Navarrete  Núcleo Milenio – Centro de Conservación Marina; Estación Costera de
Investigaciones Marinas, Pontificia Universidad Católica de Chile, Alameda
340, C.P. 6513677, Casilla 193, Correo 22, Santiago, Chile.

###### Abstract

Ecological connectivity in coastal oceanic waters is mediated by dispersion of
the early life stages of marine organisms and conditions the structure of
biological communities and the provision of ecosystem services. Integrated
management strategies aimed at ensuring long-term service provision to society
do not currently consider the importance of dispersal and larval connectivity.
A spatial optimization model is introduced to maximise the potential provision
of ecosystem services in coastal areas by accounting for the role of dispersal
and larval connectivity. The approach combines a validated coastal circulation
model that reproduces realistic patterns of larval transport along the coast,
which ultimately conditions the biological connectivity and productivity of an
area, with additional spatial layers describing potential ecosystem services.
The spatial optimization exercise was tested along the coast of Central Chile,
a highly productive area dominated by the Humboldt Current. Results show it is
unnecessary to relocate existing management areas, as increasing no-take areas
by 10% could maximise ecosystem service provision, while improving the spatial
representativeness of protected areas and minimizing social conflicts. The
location of protected areas was underrepresented in some sections of the study
domain, principally due to the restriction of the model to rocky subtidal
habitats. Future model developments should encompass the diversity of coastal
ecosystems and human activities to inform integrative spatial management.
Nevertheless, the spatial optimization model is innovative not only for its
integrated ecosystem perspective, but also because it demonstrates that it is
possible to incorporate time-varying biophysical connectivity within the
optimization problem, thereby linking the dynamics of exploited populations
produced by the spatial management regime.

Running page head: Spatial optimization of coastal eco-services

Graphical abstract.

#### Keywords:

Benthic communities, Ecological connectivity, Integrated management, Protected
areas, Aichi targets, Coastal users, Social-ecological systems, Spatial
management

## Introduction

Coastal marine ecosystems are being perturbed, fragmented, polluted, and
subjected to wholesale biodiversity loss due to a multiplicity of human
drivers, including the over-exploitation of many species of commercial
interest (Worm et al., 2006; Crain et al., 2009). Addressing this global
phenomenon requires effective implementation of sustainable management
practices that consider not only fished species, but also the multiple
ecosystem services provided by marine and coastal areas (Hutton and Leader-
Williams, 2003; Adams et al., 2004). The debate over optimal strategies to
manage and preserve marine ecosystems remains unresolved (Hilborn, 2017;
Roberts et al., 2017; Worm et al., 2006). Central to the debate is the use of
marine protected areas (MPA) as critical instruments for ecosystem
conservation and spatial management. In fact, the increased implementation of
MPAs in the last decade (Jantke et al., 2018) has been driven by international
agreements to protect 10% of Economic Exclusive Zones by 2020 (UNEP, 2010).
However, there are increasing concerns that the excessive focus on the areal
component of the 10% protection target (CBD Aichi target 11;
http://www.cbd.int/sp/targets/) is detrimental to an effective management
(Gill et al., 2017) of representative and well-connected networks of MPAs
(Visconti et al., 2019). Spatial planning must explicitly consider inherent
dependencies between marine patches under different conservation and
exploitation regimes to ensure population, genetic, community and/or ecosystem
flows or connectivity (Smith and Metaxas, 2018).

Spatial connectivity is key for the structure of marine ecosystems, as it
underlies biological as well as biophysical processes that determine ecosystem
productivity, dynamics and resilience (Gaines et al., 2010; Watson et al.,
2012; Hidalgo et al., 2017). There are different types and scales of
ecological spatial connectivity (e.g., population connectivity, genetic
connectivity, ecosystem connectivity; see Brown et al., 2016), but, in general
terms, ecological spatial connectivity refers to biological and biophysical
processes that connect areas over comparatively large spatial scales, i.e.,
those that surpass the scale of spatial management instruments and which can
be crucial for the persistence of marine populations and communities (Carr et
al., 2017). For demersal and benthic marine organisms, whose adults are
predominantly sedentary, large-scale connectivity is mediated through the
dispersal of early-life stages (i.e., eggs and larvae). Egg and larval
connectivity patterns are driven by biophysical processes characterized by
strong temporal and spatial variability (Ospina-Alvarez et al., 2018). Such
variability modulates the structure and function of benthic ecosystems, and it
is therefore an essential input for the configuration of MPA networks
(Fovargue et al., 2017). However, spatially and temporally dynamic
connectivity among marine populations is generally ignored by spatial planning
schemes, including MPA design criteria (Leslie 2005; Balbar and Metaxas 2019;
Ramesh et al. 2019; Hidalgo et al. 2019; but see, Krueck et al. 2017). Most
marine conservation planning approaches focus on "potential connectivity",
statically represented through the size and spatial placement of protected
area units (e.g., Smith and Metaxas, 2018). These approaches assume that
spatial units are self-sustained in the long-term (Pressey et al., 2007;
Magris et al., 2015). Recent models (e.g., White et al., 2014) incorporate
larval connectivity parameters, such as self-recruitment and network
centrality, into spatial optimization approaches (e.g., MARXAN and MARXAN-
Connect; Rossi et al., 2014; Dubois et al., 2016; Daigle et al., 2020) to
identify sites most relevant for meta-population persistence (Magris et al.,
2015; Krueck et al., 2017).

An additional challenge in the science of MPA network design (sensu Gaines et
al., 2010) is to determine how different configurations of protected areas
guarantee not only biological population persistence, but also maximise the
long-term provision of ecosystem services to ensure society well-being (Curtin
and Prellezo, 2010; Bennett et al., 2015). "Uses" and "activities" are
familiar concepts for marine resource managers, but, beyond direct economic
returns, the non-material benefits associated with a given use are commonly
overlooked when making decisions about protected areas (Cornu et al., 2014).
The incorporation of societies’ values into conservation proposals (Martín-
López et al., 2012; de Juan et al., 2017; Ban et al., 2019) has been slow in
terrestrial and marine habitats worldwide (Daily et al., 2009; Bennett et al.,
2015). The implementation of scientifically-based decision making that
considers feedback between ecosystems and social well-being remains a major
challenge (Congreve and Cross, 2019; Norström et al., 2020). An additional
challenge is posed by the generally small spatial scales of management
schemes, given that ecosystem services are not always demanded in the same
geographic location as the ecological processes that support them (van
Jaarsveld et al., 2005; Rodríguez et al., 2006). The biophysical connectivity
between areas is key to link ecosystems and society at regional scales (Popova
et al., 2019; Potts et al., 2014). For example, biological communities produce
larvae that can migrate long distances to recruit and replenish fished
populations (Ospina-Alvarez et al., 2015), providing resources to local
fishers. Also, carbon storage by seagrasses benefits society at a global
scale, and thus, the beneficiaries might be very distant from the ecosystem
service providers (Congreve and Cross, 2019). These processes are not embedded
in current spatial planning practices that often overlook the open nature of
marine ecosystems (Maxwell et al., 2020), e.g., by establishing artificial
frontiers in the coastal habitat continuum. Accounting for ecological
connectivity across management units at large spatial scales is therefore
crucial to ensure the provision of multiple ecosystem services in the long
term (Potts et al., 2014).

To address these multi-disciplinary challenges, we developed a spatial
optimization model that explicitly incorporates the inherent dynamism of
coastal biophysical systems (Ospina-Alvarez et al., 2018). Ecological
connectivity is estimated through coupled biophysical models of egg and larval
dispersal along a coastal system. This approach differs from most decision-
making tools by embedding larval connectivity matrices in the optimization
algorithm, rather than incorporating connectivity as a stationary property of
coastal metapopulations (e.g., MARXAN; Chollett et al., 2017). The
optimization model was built based on a case study in the central coast of
Chile, a highly productive region within the Humboldt Upwelling Ecosystem,
where wind-driven coastal upwelling supports productive industrial and small-
scale fisheries (Thiel et al., 2007; Chavez and Messié, 2009). A circulation
model is used to capture realistic dynamics of the coastal ocean in this area,
this model has previously been applied to examine dispersal and connectivity
of coastal invertebrates (Ospina-Alvarez et al., 2018; Blanco et al., 2019).
Biological data from the coastal ecosystem dominated by intertidal and
subtidal rocky habitats, the habitat of the principal commercial species for
small-scale fisheries in the region, is also available (de Juan et al., 2015).
The management regime for the small-scale fisheries follows a territorial user
rights for fisheries (TURF) spatial scheme, locally known as Areas de Manejo y
Explotación de Recursos Bentónicos (AMERBs; Fernández and Castilla, 2005;
Gelcich et al., 2012). The available data set also includes the demand for
ecosystem services by users of the coast in this area (de Juan et al., 2017).

The model prioritises marine zoning to achieve ecosystem conservation and
social benefits in areas with productive coastal ecosystems that supply a set
of services demanded by coastal users (de Juan et al., 2015). Locally, the
present work is timely because the well-established TURF system has been the
subject of heated debate in the Chilean Congress following approval of The
National Biodiversity Strategy 2017 – 2030 in 2018 (Decree 14, Ministry of
Environment, Chile). The strategy mandates the use of coastal planning
approaches that incorporate CBD recommendations for the implementation of
protected areas (Ministerio del Medio Ambiente Chile, 2018). Globally, the
proposed model offers a promising practical tool to design and implement
large-scale integrated management of connected coastal social-ecological
systems.

## Methods

### Model Domain

The model domain extends along the central coast of Chile (31.57ºS to
36.00ºS), with a spatial resolution of 2 latitudinal km (Fig.1). It
corresponds to the Eastern Boundary Upwelling Ecosystem of the south-east
Pacific (Strub et al., 1998; Thiel et al., 2007). This ecosystem is dominated
by the south-east Pacific atmospheric anticyclone, which provides increasingly
steady upwelling-favourable winds to the north, and strong but seasonally
variable winds to the south. The main surface circulation feature is the
Humboldt Current, which flows northwards, approximately 200 km off the coast.
A more coastal surface current (Chilean coastal current) also flows
predominantly to the north, and responds to upwelling-favourable wind forcing
(Aiken et al., 2008). This coastal surface current is essential for the
transport of coastal-released gametes (Ospina-Alvarez et al., 2018).

Figure 1: Study area maps showing (left) the percentage of rocky shore
substrate per unit area and (right) the percentage of managed areas (MA on
rocky shore substrate), including TURFs, Natural Monuments and Natural
Sanctuary areas.

Within the model domain there are 75 operative TURFs; 1 coastal Natural
Monument and 4 coastal Natural Sanctuary areas, covering _ca._ 31% of the
shoreline. Data describing the spatial distribution of these protected areas
(hereafter MA, e.g., Managed Areas) was obtained from the Undersecretary of
Fisheries (SUBPESCA) and the Ministry of the Environment (Ministerio del Medio
Ambiente in Spanish; Fig.1, right panel). Digitalization of the spatial data
was conducted with the software QGis v3.8. Rocky shore areas (e.g., areas
dominated by rocky vs. sandy substrate) were the target habitat in this study
as it sustains the main coastal fisheries in the region. Habitat distribution
was obtained by characterizing the habitat type (rock vs. sand) in each of the
latitudinal cells using high resolution satellite photography (Fig.1, left
panel) for the intertidal zone and projecting this onto a uniform shallow
subtidal band down to a depth of 200 m. Combined, these two data sources
describe the proportion of rocky coast covered by MA in each cell. The
remaining fraction of rocky coast in each cell was assumed to be open access
(hereafter OA; see Blanco et al., 2017, for details).

### Biological data

Shallow subtidal rocky habitats along the central Chilean coast are
characterized by kelp forests (_Lessonia trabeculata_) and a diverse
assemblage of mobile macro-invertebrates, rock fishes, birds and marine
mammals that predominantly predate or graze benthic species (Pérez-Matus et
al., 2017). Several of these rock fish and invertebrate species are subjected
to intense exploitation, as divers or coastal gatherers mostly target the
muricid gastropod, ‘loco’, exported as the Chilean abalone (p1, _Concholepas
concholepas_), key-hole limpets (p2, _Fissurella spp._) and red sea urchin
(p3, _Loxechinus albus_).

Biological surveys to estimate the net larval production (_L_) of the three
most important commercial species, loco, key-hole limpet and red sea urchin,
were conducted in the subtidal zone in 2012-2014. Four sites were selected
within the model domain to collect samples and assess the adult size and gonad
investment of key-hole limpets and red sea urchins (see Blanco et al., 2017,
for details). Data collected by the research group were used to estimate adult
size and gonad investment of Chilean abalone (n = 2,900). An additional set of
surveys obtained in seven sites within this domain, including MA, OA and no-
take areas (hereafter NT), were used to assess adult density (D; see de Juan
et al., 2015, for details). These variables (per capita gonad production and
adult density) allowed us to estimate potential fecundity (F) and, from this,
net larval production (L), which was estimated from net egg production
assuming no mortality during egg development (see eqs. 1-5 in Ospina-Alvarez
et al., 2015). Based on these studies (see Blanco et al., 2017, 2019, for
details), we estimated that net larval production in MAs was 17% less than in
NT areas, which was the baseline (Table LABEL:table:gametes). Scale factors to
estimate L in open access areas were calculated for each species from field
work and laboratory experiments (see Blanco et al., 2017). Afterwards, L per
cell was estimated based on potential fecundity per unit of rocky area for the
different fishing regimes (see Blanco et al. 2017 and eqs. from 2 to 8 in
Blanco et al. 2019). Additionally, average benthic species richness ($V_{z}$)
observed in the rocky subtidal habitat was collected by SCUBA dive surveys of
benthic communities, i.e. sub-canopy sessile benthos and mobile invertebrates,
and reef fish in the seven sampling sites, with paired MA and OA. Records of
average number of species under each fishing regimen provided estimates of an
average of 10.2 species per m2 in OA, 11.7 in MA, and 12.1 in NT (see details
on the methodology in de Juan et al. 2015).

Table 1: Scale factor for larval production (L) per species for each management regime. No-take management regime areas were considered as the baseline. |  Open Access Areas (OA) |  Managed Areas (MA) |  No-Take Areas (NT)
---|---|---|---
Loco (p1) |  0.230 |  0.830 |  1
Limpet (p2) |  0.126 |  0.830 |  1
Sea urchin (p3) |  0.350 |  0.830 |  1

### Biophysical connectivity

The ocean velocity fields were obtained from the HYbrid Coordinate Ocean Model
(HYCOM) analysis (Chassignet et al., 2007), a data assimilating forecast of
the global ocean circulation with sufficiently high resolution to accurately
reproduce mesoscale processes that dominate ocean variability. Then, a
Spatially Explicit Individual Based Model (SEIBM) was used to simulate the
dispersal of early life stages of the three commercial species. The SEIBM was
coupled with the 3D hydrodynamic model using a customized version of the open
source modelling tool ICHTHYOP (Lett et al., 2008). The code was used to
simulate trajectories of early life stages from velocity fields in a 3D
hydrodynamic model, using Pelagic Larval Durations (PLD) of the three model
species: 90 days for the Chilean abalone, 20 days for the red-sea urchin and
10 days for key-hole limpets. During this development time, larvae are
subjected to advective and diffusive processes that condition their
alongshore, cross-shore and vertical movement in the water column and,
consequently, the total "larval waste" and successful onshore recruitment at
the end of the development period. In consequence, larvae leaving the studied
domain through the northern, southern and western open boundaries and larvae
advected into oceanic waters are considered to be lost for recruitment. We
considered that larvae display Diel Vertical Migration (DVM) behaviour,
previously shown to effect recruitment rates onshore but which does not
significantly alter spatial or temporal patterns of connectivity (Ospina-
Alvarez et al., 2018). A connectivity score that incorporates seasonality was
used to capture observed reproductive peaks of the model species: summer and
winter for _Concholepas_ , and spring and summer for key-hole limpets and the
red sea urchin. The connectivity score represents the standardized recruitment
number of particles transported from spawning to recruitment locations
weighted by demographic and habitat data. The spawning locations were assumed
to be dependent only on available habitat and were therefore identical for all
species. The optimal number of particles to be released in order to reach
stable patterns was established at 20,000 per day, every five days (i.e.,
120,000 particles per month). For a detailed description of the methodology
see Ospina-Alvarez et al. (2018).

A "connectivity matrix" (C; as defined by Cowen and Sponaugle 2009), per
dispersal event was generated. Potential connectivity was defined as the
probability of larval transport from a spawning site _j_ (columns in C) to a
destination location _i_ (rows in C), after completing development. Thus, the
diagonal of C is the probability of local larval retention, and the sub-
diagonals indicate the probability of dispersal between sites. The C matrix
contains all the necessary information for metapopulation dispersal dynamics
and time-varying processes (Aiken and Navarrete, 2014): the character of the
ocean dictates that it is typically time-varying and asymmetrical. Realized
larval connectivity corresponds to larval supply and it can be defined as the
effective number of larvae that travel from _j_ to _i_ and the corresponding
matrix ($H_{jip}$) can be estimated for each species (_p_) using potential
connectivity matrices weighted by relevant biological and environmental
information (Watson et al., 2010). Specifically, habitat availability
(percentage of rocky and sandy coastline), management regime (percentage of
MA, NT and OA) and adult abundance and size (as a factor determining egg and
larval production) per cell were used to weight the potential matrix. All
these factors were considered invariant over the timescales of dispersal.

In consequence, the model was structured into six main components:

  * •

environmental envelope: proportion of rocky coast and prevailing management
regime $z\;\in\;\left(OA,\;MA,\;NT\right)\ $.

  * •

source sites (i.e., _j_ , spawning sites).

  * •

destination sites (_i_ , settlement sites).

  * •

net larval production (L) as a function of density, size and fecundity of
adult individuals.

  * •

potential larval connectivity, C: probability that larvae spawned at site, _j_
, will end up at site _i_.

  * •

realized larval supply (H): C matrix multiplied by the number of larvae
released in cell _j_ which is a function of the larval per capita production,
the environmental envelope of the source cell (binary) and the abundance of
adults in that cell. The management regimen does not play a role as
microscopic larvae are transported throughout the study site (Fig.2).

Figure 2: Larval connectivity networks for key-hole limpet (p2), red sea
urchin (p3) and _C. concholepas_ (p1). The connection between two locations is
represented as curved lines (arcs) in a yellow-red colour scale that
represents larval connectivity probability. Recruitment intensity is
represented by dots (green-blue colour scale). The larval connectivity
probability is extracted from C matrices and the recruitment index represents
the standardized number of recruited larvae at each location. The figure is
based on patterns and results described by Ospina-Alvarez et al. (2018).

### Quantification of ecosystem services

The prioritisation of ecosystem attributes in the study area was assessed
through face-to-face interviews with the principal end-users (de Juan et al.,
2017). Over 900 questionnaires were conducted in six coastal locations within
the model domain, where small-scale fisheries and national tourism are the
principal activities; therefore, these questionnaires targeted local fishers,
residents and tourists. The survey was designed to elicit user groups’
prioritisation of ecosystem attributes. Despite socioeconomic variability
among the study locations, ecosystem valuation by users, and among users’
groups, was highly homogeneous. Several attributes were relevant to all users,
like scenic beauty and clean waters and beaches; however, fishers also
prioritised fishing resources and biodiversity, while tourists and residents
prioritised intangible values including the peace and relaxing atmosphere of
the coast (see more details on the approach in de Juan et al., 2017).

The potential of coastal areas to provide these services was estimated from
bio-physical data in the area. The highly valued "scenic beauty" of the
coastline was identified by users as the absence of urban development on rocky
coasts (de Juan et al., 2017). This attribute can also be indirectly related
with clean water, and peace and relaxation, the other intangible services
valued by users. The proportion of the coast within the model domain that has
not been modified by urban development was identified through Google Earth
images. This information was overlapped with the characterization of the coast
as rocky substrate and provided the proportion of each cell characterized by
the attribute scenic beauty. The biological diversity and fishing resources,
mainly prioritised by fishers, are also are included in our model.

The connectivity matrices C, described in the previous section, define which
areas make disproportionately large contributions to realized larval supply to
other areas along the coast. Sites that excel for their net larval production
ultimately supply adult populations in the region and are key to maintaining
populations of commercial species in the area, as long as they themselves
receive enough individuals to maintain the local population. Therefore,
biological productivity in the region is assumed to be well-correlated with
net larval production at a given site, since this variable integrates adult
abundance and size. Based on its link with biological productivity, net larval
production was used as a proxy for the potential of a site to secure ”food
production” (i.e., cath of fishing resources) as an ecosystem service.

The average number of benthic species ($V_{z}$) estimated for the coastal
areas under the different protection regimens was used as an additional input
variable in the optimization model, when the target was to maximise the
biological diversity of a site (see next section for more details) as the
baseline for the provision of numerous ecosystem services, including
recreational value (as identified in de Juan et al., 2017).

### Spatial optimization model

We developed a spatial optimization model to identify the optimal
configuration of marine zoning (e.g. locations of MA, NT and OA in the study
area) that maximises the potential of coastal areas to meet integrated
ecosystem management targets. Our modelling framework builds on previous
approaches (e.g., Davis et al., 2015; Watts et al., 2009), but uses positive
decision variables with linear programming versus binary decision variables
(using, e.g. mixed integer programming). All models are solved in GAMS
(General Algebraic Modelling System, v.23.5.2), available at
https://www.gams.com/, using the CPLEX solver (CPLEX, 2009). We considered
three targets, which aim to maximise the potential provision of different sets
of ecosystem services. In target 1 (T1), we maximise net larval production (as
a proxy for food production). In target 2 (T2), we maximise net larval
production and biodiversity. In target 3 (T3), we maximise net larval
production, biodiversity, and scenic beauty. Spatial priorities consistent
with target 1 require identifying the best spawning j sites in terms of net
larval production, weighted by the proportion of larvae successfully arriving
to any site _i_ , for each commercial species. Biodiversity (T2 and T3) is
represented through benthic species richness. To increase biodiversity in the
study region, spatial zoning should prioritise zones that support higher
levels of benthic species richness, e.g., NT. Finally, scenic beauty (T3)
values are increased by locating NT areas in parts of the coast with high
scenic beauty values, e.g., cells where > 80 % of total area is rocky habitat
and there is no urban development.

For each target, we assessed a range of scenarios based on existing or
potential management regimes in central Chile. Scenario 1 (S1) is our baseline
scenario, and describes the existing management regime in central Chile. This
scenario provides a baseline against which we explored the potential for
alternative management scenarios to increase (or decrease) the potential
provision of ecosystem services to society. Scenarios 2 and 5 describe
pessimistic and optimistic (S2 and S5 respectively) outcomes for potential
ecosystem service provision in the study region. Scenario 3 is based on
realistic projections of management changes that could take place in the
region. Finally, scenario 4 incorporates the objectives of the CBD targets by
considering the assignation of 10% of the coast to no-take areas based on a
spatial configuration that aims to maximise biological connectivity
constrained by holding an integral approach that explicitly considers multiple
ecosystem attributes.

  * •

Scenario 1 (S1): Maintains existing management conditions.

  * •

Scenario 2 (S2): Final allocation is constrained to 100 % OA areas. This is a
worst-case scenario from the point of view of conservation and sustainability.

  * •

Scenario 3 (S3): Model is constrained by existing ratio of MA and OA areas
(30:70) but is free to allocate these areas in space.

  * •

Scenario 4 (S4): Model solution is constrained to include existing MA, but
freely allocates 10% of the rocky study area to NT (from existing OA areas).

  * •

Scenario 5 (S5): Final allocation is constrained to 100 % NT areas. i.e. there
is loss of an ecosystem service (fishery), but it achieves the maximum
possible conservation.

The characteristics of each scenario constrain the integrated ecosystem
management target that can be achieved. Net larval production can be maximised
(or assessed) across all five scenarios, therefore we assess the optimal
configuration of the study region to meet target 1 for each scenario. However,
for targets 2 and 3 we consider only a subset of scenarios. To achieve target
2, e.g., maximise net larval production and biodiversity, we assess scenarios
3-5. Scenarios 1-2 are not assessed for this target because they contain no NT
zones, therefore outcomes would repeat the outcomes of target 1. Target 3
maximises net larval production, biodiversity and scenic beauty. As scenic
beauty is maximised by allocating a cell with rocky habitat, and >80 % of its
area free from urban development, to NT, we limit analysis for this target to
scenario 4.

To identify the optimal spatial configuration of the study area that meets
each integrated ecosystem management target, we define three objective
functions. The constraints applied to each function are specific to the five
scenarios. We use linear programming with positive decision variables
($X_{iz}$) to identify the optimal allocation of each cell _i_ to each
management zone _z_ (MA, NT, OA). Decision variable $X_{iz}$ will take a
positive value if cell _i_ is allocated to management regime _z_ , and 0
otherwise. The following constraint (eq. 1) ensures that each cell is fully
allocated across the three zones, but not over-allocated (e.g. > 1).

$\sum_{iz}{X_{iz}\;=\;1}$ (1)

The objective function for target 1 (e.g. maximise net larval production)
incorporates realized larval supply and the probability that a larvae will
settle in a given cell. As defined previously, realized larval supply
($H_{jip}$), describes potential larval connectivity, the number of larvae
released from a destination cell (_j_), and other qualities of the destination
cell, including the environmental envelope and abundance of adults of a given
species (_p_). The probability that a larvae will settle in a given cell is
determined by settlement probability $\lambda_{ipz}$, which is specific for
each species _p_ and each management regime _z_ , and the environmental
envelope of the destination site, i.e., the percentage of available rocky
habitat in each destination cell _i_. We maximise the following objective
function (Eq. 2) to achieve target 1:

$\max\sum_{ipz}{X_{iz}\;H_{jip}\;\lambda_{ipz}}$ (2)

In scenarios 3 and 4, this objective function is subject to an additional
constraint, which ensures that 30 % of the rocky area ($R_{i}$) in the study
region is allocated to the MA zone (S3), and 10 % of rocky area is allocated
to the NT zone (S4):

$\sum_{iz}{R_{i}}\;\times\;M_{z}\;=\;\sum_{iz}{X_{iz}}R_{i}$ (3)

Specifically, this constraint (Eq. 3) specifies that the rocky area in each
cell allocated to a specific management regime, must be equal to $M_{z}$,
where $M_{MA}$ = 0.3 for scenario 3, and $M_{NT}$ = 0.1 for scenario 4.

The objective function for target 2 builds on the objective function of target
1, including all general (Eq. 1), and scenario-specific constraints (Eq. 3),
but includes a term to describe the diversity value $V_{z}$ of each management
regime _z_ (eq. 4).

$\max\sum_{ipz}{X_{iz}\;H_{jip}\;\lambda_{ipz}\;V_{z}}$ (4)

The objective function for target 3 is the same as target 2 but subject to an
additional constraint: that cells allocated to the NT zone must have at least
80% of their rocky area free from human construction.

We normalised the following $x\;\in\;\left(H_{jip},\;V_{z}\right)$ ecosystem
service values between 0 and 1 following the equations below (eqs. 5-7), so
that realized larval production and diversity values were comparable and
contributed equally to the optimization value.

$a\;=\;min(x)$ (5)

$b\;=\;max(x)$ (6)

$x_{n}\;=\;\frac{\left(x-a\right)}{\left(b-a\right)}$ (7)

## Results

Net larval production slightly increased under the existing management
scenario in central Chile (S1), compared to the worst-case scenario with 100%
of the shore allocated to OA areas (S2; Table 2). However, by allowing the
optimization model to freely allocate the spatial location of the existing MA
and OA areas (S3), i.e. without altering total area (and therefore stock)
proportions, the objective function (17.49) increased relative to S1 and S2
(Table 2 or Fig.3). A similar objective function value (17.94) was obtained
when an additional 10% of open access area was re-allocated to no-take area
(S4), without altering the existing geographic distribution of MA. These
scenarios (S3 and S4) represent a 3- and 5-fold improvement over the current
and worst-case scenarios, respectively. However, under these conditions, the
management areas (MA) and the no-take areas (NT) are highly overlapped and
concentrated in the central and northern sections of the model domain, which
may be undesirable (Fig.4, panel T1S3 and T1S4). As expected, allocating the
entire shore to no-take areas (S5) produced the highest net larval production
(objective function: 39.08). It is important to emphasize that net larval
production is being optimized by embedding the connectivity matrices in the
optimization algorithm. Therefore, a site with low connectivity or with low
self-recruitment will not get a high priority in the optimization exercise,
even if that site had a high net larval production.

Table 2: Targets, scenarios and optimization values. Net larval production
values are described in Eqs. 2 and 4), species diversity optimization value
corresponds to $V_{z}$ in (Eq. 4). The objective function is the resulting
output from the optimization equation and constraints applied to each target-
scenario combination. Consider that diversity ($V_{z}$) for T2 and T3 is an
average value for each type of management area _z_ (MA, NT, OA) and that
scenic beauty is a constraint only affecting the allocation in NT for T3. In
consequence, some target-scenario combinations are not possible or have the
same value as scenarios evaluated under T1.

Target number |  |  | Optimization value
---|---|---|---
and description | Scenario number | Scenario description |  | Larval
---
connectivity
| Species
---
diversity
| Objective
---
function
1 Net larval production (NLP) | 1 | Existing | 5.96 | - | 5.96
2 | 100% allocated to OA | 3.45 | - | 3.45
3 |  | Existing MA proportion but optimal
---
allocation [MA: 0.3, OA: 0.7]
17.49 | - | 17.49
| 4 |  | Existing MA allocation, but additional
---
10% of study area changed from
OA to NT zones [MA: 0.3, OA: 0.6, NT: 0.1]
17.94 | - | 17.94
| 5 | 100% allocated to NT | 39.08 | - | 39.08
2 Net larval production + Biodiversity (NLP + Biodiversity) | 3 |  | Existing MA proportion but optimal
---
allocation [MA: 0.3, OA: 0.7]
17.49 | 51.21 | 16.95
4 |  | Existing MA allocation, but additional
---
10% of study area changed from
OA to NT zones [MA: 0.3, OA: 0.6, NT: 0.1]
17.74 | 51.63 | 17.65
| 5 | 100% allocated to NT | 39.08 | 55.99 | 39.08
| 3
---
Net larval production +
Biodiversity + Scenic beauty
(NLP + Biodiversity + Scenic beauty)
4 |  | Existing MA allocation, but additional
---
10% of study area into
NT (area taken from OA) [MA: 0.3, OA: 0.6, NT: 0.1]
15.43 | 51.63 | 15.33

Figure 3: Optimization results for all possible combinations of targets and
scenarios. Notice that some target-scenario combinations are not possible or
have the same value as scenarios evaluated under Net larval production (NLP,
T1) because for NLP + Biodiversity (T2) and NLP + Biodiversity + Scenic beauty
(T3) species diversity ($V_{z}$) is an average value for each type of
management area _z_ (MA, NT, OA) and scenic beauty is a constraint only
affecting the allocation in NT for T3. Figure 4: Scenarios to allocate
management or/and no take areas from the optimization model. Net larval
production – existing MA (T1S1); Net larval production – Existing MA
proportion but optimal allocation (T1S3); and Net larval production – Existing
MA allocation, but additional 10% of study area changed from OA to NT zones
(T1S4). The legends for % MA and geographic features and the scale are shared
by all panels.

In scenarios 3 and 4, the addition of biological diversity to net larval
production (T2) obtained a similar objective function value as target 1 (i.e.,
only considering net larval production) (Table 2). This result suggests that
the simultaneous maximisation of diversity and net larval production was not
detrimental to net larval production in this ecosystem (Fig.3), probably
because MA are located in productive areas with high biological diversity. In
consequence, scenarios 3 and 4 for target 1 and target 2 are very similar and
most NT are allocated to cells with a proportion of MA, resulting in numerous
small NT (Fig.5, panel T2S4).

Figure 5: Scenarios to allocate management or/and no take areas from the
optimization model. Net larval production + Biodiversity – Existing MA
proportion but optimal allocation (T2S3); Net larval production + Biodiversity
– Existing MA allocation, but additional 10% of study area changed from OA to
NT zones (T2S4); and Net larval production + Biodiversity + Scenic beauty –
Existing MA allocation, but additional 10% of study area changed from OA to NT
zones (T3S4). The legends for % MA and geographic features and the scale are
shared by all panels; the legend for % NT is shared by T2S4 and T3S4 panels.

When we incorporate all the variables prioritised by end-users, by the
inclusion of scenic beauty into conservation targets (target 3), and we
increase 10% of no-take area, there was a slight decrease in optimization
value, with respect to target 2 (T2S4), from 17.65 to 15.33 (Table2). This is
a consequence of constraining the best areas for biodiversity and larval
connectivity to the presence of "scenic beauty". In this scenario there was an
improved spatial representativity of NT, with less overlap of NT and MA as
observed in T2S4, as cells with 100% of OA are allocated to NT (Fig.5, panel
T3S4). Note that all these scenarios provide larger optimization functions
than the existing scenario.

## Discussion

The spatial optimization model introduced in this work addresses the challenge
of achieving a diversity of management objectives, ranging from conservation
of coastal biodiversity, to the maximisation of net larval production and
larval connectivity of key commercial species for local fisheries, and to
societal preferences for the scenic beauty of the coast. The application of
this novel integrated approach to the central coast of Chile is promising, and
highlights that an additional 10% protection of the coast could increase
5-fold the potential for ecosystem services’ provision. The spatial
optimization targets were founded on the basis that the fisheries management
regimen (MA, NT, OA) will change the population density and size of commercial
species, and the biological diversity. These changes affect the realized
larval supply (proxy of biological productivity) along the coast mediated by
biophysical connectivity. The different targets provided similar outputs,
likely driven by correlation between variables, e.g., biodiversity and larval
connectivity or biodiversity and scenic beauty provided by non-urbanized
areas, which is a consequence of exploring the value of coastal areas already
shaped by human intervention. Despite limited variability in the overall
output (i.e., the objective function), the simultaneous achievement of
multiple conservation objectives (target 3) implied a different geographical
distribution of protected areas compared to the current scenario or to target
1 (when a single objective was considered), with large sections of OA areas
assigned to NT, resulting in lower number of NT. This solution might be more
feasible than small NT distributed over larger number of cells, which might
imply high costs in surveillance. Regardless of the target, an optimal
scenario implied maintaining the current location of MA but the addition of
10% of the coast to NT. Thus, along our modelled coastline, it would be
possible to achieve integrated conservation objectives and reduce the
perceived negative effects of MPA implementation (Ban et al., 2019) at no
additional cost for fishers in terms of the removal of users’ rights in
existing management areas. Still, fishers opposition might be detected for the
closing of large OA areas (as in target 3), whereas smaller NT adjacent to MA
(as in targets 1 and 2) could be supported by local fishers as part of their
fisheries management plan. It is important to note that a case-by-case
approach must be followed to explore alternatives of no-take area positioning
that minimizes other perceived negative impacts on end-users (Ban et al.,
2019).

The optimization exercise constrained by a target that incorporated scenic
beauty along the coast (i.e., non-urbanized areas as a proxy for scenery,
clean beaches and peaceful environments, as prioritised by end-users in the
area, (de Juan et al., 2017)) illustrates the multi-dimensional nature of
coastal ecosystem services and the feasibility of simultaneously including a
set of very different ecosystem services in management decisions (e.g.,
provision, regulating and cultural services, Millennium Ecosystem Assessment,
2005). Despite the simultaneous inclusion of three conservation targets
causing a decrease in the objective value, it further improved the spatial
representation (spacing) of protected areas (NT and MA) along the coast. It
must be noted that the spatially uneven distribution of restricted areas
across scenarios is a limitation of the model, as it prioritises regions where
the rocky habitat for the commercial species is dominant. Previous studies
have criticized the lack of analysis assessing the extent to which the
representativity of many sensitive habitats has been achieved within global
MPAs (Fischer et al., 2019; Jantke et al., 2018). Our current approach aims to
improve habitat connectivity and representativity in a network of fishery
restricted areas that are connected through the strong hydrodynamism of the
Humboldt current (Ospina-Alvarez et al., 2018). However, the model was
designed for the rocky coast, meaning that the central region was consistently
prioritised as a conservation area, as the southern section of the model
domain is characterized by greater proportions of sandy coast. This result
emphasizes that habitat representativeness is an important criterion for
reserve network design. Further model development should incorporate key
species and diversity values from the sandy subtidal areas in the region, so
habitat diversity is considered. This model extension would also allow
incorporation of additional societal preferences in the central coast of
Chile, linked with the diversity of coastal habitats, e.g., clam fisheries,
beach tourism, surf, resulting in a truly integrative spatial optimization
approach.

The positioning of management areas was highly consistent across scenarios and
there is probably an historical explanation for this, as fishers tend to
select the most productive sites for the allocation of management areas.
Moreover, these productive sites, with strong larval connectivity, are
probably also high diversity sites as other species would also match the
connectivity patterns of the key species. Therefore, some management areas
currently located in the best places to provide ecosystem services could be
important ancillary conservation instruments if they were well enforced. The
level of enforcement is key, as previous studies in central Chile observed
that densities of economically important macro-invertebrates and reef fish are
not significantly different between highly enforced management areas and no-
take MPAs (Gelcich et al., 2012). These observations are critical as from the
twenty official MPAs in Chile, only five are effectively managed (Mora et al.,
2006; Petit et al., 2018). During the last decade, Chile has made some
advances in the establishment of protected areas. However, a pending task is
to consolidate an MPA network in the coastal area, particularly in the central
region where human-uses are concentrated and no-take areas would probably face
fishers opposition (Suman et al., 1999; Gelcich et al., 2009). In this
context, the current exercise aimed to develop a tool that could assist
prioritisation of management efforts for the location of protected areas and
the identification of management areas that could incorporate a no-take
section and/or an area to be managed under special regimes to increase
enforcement level. Priority sites for conservation frequently overlapped with
highly populated areas, which include large commercial harbours and polluted
areas associated with hydrothermal plants (e.g., central prioritised cells in
the model overlap with the urban area of Valparaiso). This overlap implies an
unfeasible establishment of protected areas and highlights the need to
consider all critical components of coastal ecosystems in the spatial
optimization and increase spatial resolution of the information and models.
With the incorporation of scenic beauty in the optimization exercise, we
largely overcome this issue as, by prioritising sites that hold higher scenic
beauty, it indirectly selects less urbanized sites. These local issues, that
complicate spatial planning approaches, are common to many countries that lack
a legal framework for the spatial management of the coastal and marine
continuum. In these scenarios, the identification of priority sites relying on
multiple and interacting data sources is essential for decision making.

By incorporating ecological connectivity in the spatial prioritisation of
marine resources, we show that the effects of management efforts can propagate
beyond the region through biophysical processes that support the provision of
ecosystem services by coastal areas. Continuous improvement of physical and
biological information is desirable and necessary before any geographically-
specific implementation of model results. For instance, higher-resolution
information on biological community structure, population size distribution
and diversity, including habitat diversity (e.g., subtidal sands or kelp
forests), is necessary to capture small-scale biological spatial variability.
Currently, the model considers average biological metrics on "rocky shores"
per management regime, shaped by the dynamism of a mid-resolution biophysical
model that creates patterns of connectivity. After the connectivity modelling
results, sensitivity analyses are also necessary to evaluate the different
causes of larval losses and compare them with the dynamics of larval
recruitment under natural conditions. In this work, we found that losses of
larvae to the northern and southern boundaries and offshore western were
moderate (57.6% SD ± 21.3). However, most of them were related to oceanic
advection, while losses to the North were low and occurred mainly in the
austral summer and especially for long PLDs. On the other hand, losses due to
beaching were moderate (39.3% SD ± 27.3). Total mortality, the sum of all
sources of losses, was 81.7% SD ± 19.2. The low percentages of losses towards
the northern and southern boundaries suggest that the high density of
recruitment in the central zone of the study area is not an artefact of the
larval transport model (i.e., SEIBM) but it is related to the
representativeness of the rocky habitat. The central zone of the study area is
the zone with the highest percentage of rocky habitat per cell. In our model,
a condition for larval recruitment is the existence of rocky habitat,
simulating the selection of settlement habitat by organisms living in natural
conditions. An additional issue that needs to be resolved is the incorporation
of temporal variability in the optimization exercise. The operational model on
which implementation is based should deal with temporal variability in the
hydrodynamic connectivity at scales relevant for management.

Our model was restricted to the rocky subtidal and the values linked to this
coast, which can exclude social-economic important activities such as beach
tourism or clam fisheries. Future model developments should encompass the
diversity of coastal ecosystems and human activities to inform a truly
integrative spatial management. To achieve this, research efforts in the area
should focus on acquiring biophysical data at continuous spatial scales as,
for example, another model limitation was biodiversity values included as
static variables linked to the presence or absence of fisheries restriction.
We urge the compilation of spatially-explicit databases to obtain more precise
estimates of the critical processes that sustain social-ecological systems in
coastal areas. Through sensitivity analysis, the framework proposed here can
be used to identify the type of information and the geographic location where
new data would most contribute to improve model results. The broad adoption of
ecosystem service science in coastal management is limited by the paucity of
data (Cárcamo et al., 2013; Saunders et al., 2016), therefore, models such as
ours are designed to motivate the acquisition of critical information. These
models should work in an adaptive way by gradually improving with the
incorporation of new knowledge. The simultaneous achievement of multi-
dimensional objectives is a global and critical issue for the spatial planning
of human activities on the coasts, and our modelling approach offers a path in
this direction. A wider application of the model includes the creation of
short-, medium- and long-term projections for policy makers to promote the
sustainable use of coastal ecosystems more effectively at a regional scale.
The model could easily accommodate IPCC projections of sea-surface temperature
and ocean currents to forecast larval growth and dispersal pathways coupled
with shared social-economic scenarios to analyse feedback between climate
change and social well-being (O’Neill et al., 2014). Moreover, the
visualisation of optimal configurations of management units in maps could also
assist an effective communication with relevant stakeholders.

## Author Contributions

AO, SdJ and KD conceived the ideas and designed methodology; AO, SdJ, MF and
CG contributed to data; AO, KD and SN performed model simulations and formal
analysis; AO and SdJ led the writing of the manuscript. All authors
contributed critically to the drafts and gave final approval for publication.

## Acknowledgements

This work was supported by the Millennium Scientific Initiative ICM from
Ministry of Economy, Development and Tourism of Chile to the Center for Marine
Conservation [CCM RC 1300004]; AO was supported by the National Fund for
Scientific and Technological Development, FONDECYT with a post-doctoral grant
[grant number: 3150425] and by H2020 Marie Skłodowska-Curie Actions [746361].
KD Was supported by the Australian Research Council Centre of Excellence for
Environmental Decisions [CE11001000104], funded by the Australian Government.
Additional funding was provided by Fondecyt # 1160289 to SAN. Funding for the
development of HYCOM has been provided by the National Ocean Partnership
Program and the Office of Naval Research. Data assimilative products using
HYCOM are funded by the U.S. Navy. Computer time was made available by the DoD
High Performance Computing Modernization Program. The output is publicly
available at http://hycom.org/

## References

  * Adams et al. (2004) Adams, W.M., Aveling, R., Brockington, D., Dickson, B., Elliot, J., Hutton, J., Roe, D., Vira, B., Wolmer, W., 2004.  Biodiversity conservation and the eradication of poverty.  Science 306, 1146–1149.
  * Aiken et al. (2008) Aiken, C.M., Castillo, M.I., Navarrete, S.A., 2008.  A simulation of the Chilean Coastal Current and associated topographic upwelling near Valparaíso, Chile.  Continental Shelf Research 28, 2371–2381.
  * Aiken and Navarrete (2014) Aiken, C.M., Navarrete, S.A., 2014\.  Coexistence of competitors in marine metacommunities: environmental variability, edge effects, and the dispersal niche.  Ecology 95, 2289–2302.
  * Balbar and Metaxas (2019) Balbar, A.C., Metaxas, A., 2019\.  The current application of ecological connectivity in the design of marine protected areas.  Global Ecology and Conservation 17, e00569.
  * Ban et al. (2019) Ban, N.C., Gurney, G.G., Marshall, N.A., Whitney, C.K., Mills, M., Gelcich, S., Bennett, N.J., Meehan, M.C., Butler, C., Ban, S., Tran, T.C., Cox, M.E., Breslow, S.J., 2019.  Well-being outcomes of marine protected areas.  Nature Sustainability 2, 524–532.
  * Bennett et al. (2015) Bennett, E.M., Cramer, W., Begossi, A., Cundill, G., Díaz, S., Egoh, B.N., Geijzendorffer, I.R., Krug, C.B., Lavorel, S., Lazos, E., Lebel, L., Martín-López, B., Meyfroidt, P., Mooney, H.A., Nel, J.L., Pascual, U., Payet, K., Harguindeguy, N.P., Peterson, G.D., Prieur-Richard, A.H., Reyers, B., Roebeling, P., Seppelt, R., Solan, M., Tschakert, P., Tscharntke, T., Turner II, B.L., Verburg, P.H., Viglizzo, E.F., White, P.C., Woodward, G., 2015.  Linking biodiversity, ecosystem services, and human well-being: three challenges for designing research for sustainability.  Current Opinion in Environmental Sustainability 14, 76–85.
  * Blanco et al. (2017) Blanco, M., Ospina-Alvarez, A., González, C., Fernández, M., 2017\.  Egg production patterns of two invertebrate species in rocky subtidal areas under different fishing regimes along the coast of central Chile.  PLoS One 12, e0176758.
  * Blanco et al. (2019) Blanco, M., Ospina-Alvarez, A., Navarrete, S.A., Fernández, M., 2019\.  Influence of larval traits on dispersal and connectivity patterns of two exploited marine invertebrates in central Chile.  Marine Ecology Progress Series 612, 43–64.
  * Brown et al. (2016) Brown, C.J., Harborne, A.R., Paris, C.B., Mumby, P.J., 2016\.  Uniting paradigms of connectivity in marine ecology.  Ecology 97, 2447–2457.
  * Cárcamo et al. (2013) Cárcamo, P.F., Garay-Flühmann, R., Gaymer, C.F., 2013.  Opportunities and constraints of the institutional framework for the implementation of ecosystem-based management: The case of the Chilean coast.  Ocean & Coastal Management 84, 193–203.
  * Carr et al. (2017) Carr, M.H., Robinson, S.P., Wahle, C., Davis, G., Kroll, S., Murray, S., Schumacker, E.J., Williams, M., 2017\.  The central importance of ecological spatial connectivity to effective coastal marine protected areas and to meeting the challenges of climate change in the marine environment.  Aquatic Conservation-Marine And Freshwater Ecosystems 27, 6–29.
  * Chassignet et al. (2007) Chassignet, E.P., Hurlburt, H.E., Smedstad, O.M., Halliwell, G.R., Hogan, P.J., Wallcraft, A.J., Baraille, R., Bleck, R., 2007\.  The HYCOM (HYbrid Coordinate Ocean Model) data assimilative system.  Journal of Marine Systems 65, 60–83.
  * Chavez and Messié (2009) Chavez, F., Messié, M., 2009\.  A comparison of Eastern Boundary Upwelling Ecosystems.  Progress in Oceanography 83, 80–96.
  * Chollett et al. (2017) Chollett, I., Garavelli, L., O’Farrell, S., Cherubin, L., Matthews, T.R., Mumby, P.J., Box, S.J., 2017.  A Genuine Win-Win: Resolving the “Conserve or Catch” Conflict in Marine Reserve Network Design.  Conservation Letters 10, 555–563.
  * Congreve and Cross (2019) Congreve, A., Cross, I.D., 2019\.  Integrating ecosystem services into environmental decision-making.  Journal of Applied Ecology 56, 494–499.
  * Cornu et al. (2014) Cornu, E.L., Kittinger, J.N., Koehn, J.Z., Finkbeiner, E.M., Crowder, L.B., 2014.  Current practice and future prospects for social data in coastal and ocean planning.  Conservation Biology 28, 902–911.
  * Cowen and Sponaugle (2009) Cowen, R.K., Sponaugle, S., 2009\.  Larval dispersal and marine population connectivity.  Annual review of marine science 1, 443–466.
  * CPLEX (2009) CPLEX, I.I., 2009.  V12. 1: User’s Manual for CPLEX.  International Business Machines Corporation 46, 157.
  * Crain et al. (2009) Crain, C.M., Halpern, B.S., Beck, M.W., Kappel, C.V., 2009\.  Understanding and Managing Human Threats to the Coastal Marine Environment.  Annals of the New York Academy of Sciences 1162, 39–62.
  * Curtin and Prellezo (2010) Curtin, R., Prellezo, R., 2010\.  Understanding marine ecosystem based management: A literature review.  Marine Policy 34, 821–830.
  * Daigle et al. (2020) Daigle, R.M., Metaxas, A., Balbar, A.C., McGowan, J., Treml, E.A., Kuempel, C.D., Possingham, H.P., Beger, M., 2020\.  Operationalizing ecological connectivity in spatial conservation planning with Marxan Connect.  Methods in Ecology and Evolution , 185\.
  * Daily et al. (2009) Daily, G.C., Polasky, S., Goldstein, J., Kareiva, P.M., Mooney, H.A., Pejchar, L., Ricketts, T.H., Salzman, J., Shallenberger, R., 2009.  Ecosystem services in decision making: time to deliver.  Frontiers in Ecology and the Environment 7, 21–28.
  * Davis et al. (2015) Davis, K., Kragt, M., Gelcich, S., Schilizzi, S., Pannell, D., 2015.  Accounting for enforcement costs in the spatial allocation of marine zones.  Conservation Biology 29, 226–237.
  * Dubois et al. (2016) Dubois, M., Rossi, V., Giacomi, E.S., Haond, S.A., López, C., García, E.H., 2016\.  Linking basin-scale connectivity, oceanography and population dynamics for the conservation and management of marine ecosystems.  Global Ecology and Biogeography 25, 503–515.
  * Fernández and Castilla (2005) Fernández, M., Castilla, J.C., 2005\.  Marine Conservation in Chile: Historical Perspective, Lessons, and Challenges.  Conservation Biology 19, 1752–1762.
  * Fischer et al. (2019) Fischer, A., Bhakta, D., Macmillan-Lawler, M., Harris, P., 2019\.  Existing global marine protected area network is not representative or comprehensive measured against seafloor geomorphic features and benthic habitats.  Ocean & Coastal Management 167, 176–187.
  * Fovargue et al. (2017) Fovargue, R., Bode, M., Armsworth, P.R., 2017.  Size and spacing rules can balance conservation and fishery management objectives for marine protected areas.  Journal of Applied Ecology 55, 1050–1059.
  * Gaines et al. (2010) Gaines, S.D., White, C., Carr, M.H., Palumbi, S.R., 2010\.  Designing marine reserve networks for both conservation and fisheries management.  Proceedings of the National Academy of Sciences of the United States of America 107, 18286–18293.
  * Gelcich et al. (2012) Gelcich, S., Fernández, M., Godoy, N., Canepa, A., Prado, L., Castilla, J.C., 2012\.  Territorial User Rights for Fisheries as Ancillary Instruments for Marine Coastal Conservation in Chile.  Conservation Biology 26, 1005–1015.
  * Gelcich et al. (2009) Gelcich, S., Godoy, N., Castilla, J.C., 2009.  Artisanal fishers’ perceptions regarding coastal co-management policies in Chile and their potentials to scale-up marine biodiversity conservation.  Ocean & Coastal Management 52, 424–432.
  * Gill et al. (2017) Gill, D.A., Mascia, M.B., Ahmadia, G.N., Glew, L., Lester, S.E., Barnes, M., Craigie, I., Darling, E.S., Free, C.M., Geldmann, J., Holst, S., Jensen, O.P., White, A.T., Basurto, X., Coad, L., Gates, R.D., Guannel, G., Mumby, P.J., Thomas, H., Whitmee, S., Woodley, S., Fox, H.E., 2017\.  Capacity shortfalls hinder the performance of marine protected areas globally.  Nature 543, 665–669.
  * Hidalgo et al. (2017) Hidalgo, M., Kaplan, D.M., Kerr, L.A., Watson, J.R., Paris, C.B., Browman, H.I., 2017\.  Advancing the link between ocean connectivity, ecological function and management challenges.  ICES Journal of Marine Science 74, 1702–1707.
  * Hidalgo et al. (2019) Hidalgo, M., Rossi, V., Monroy, P., Ser-Giacomi, E., Hernández-García, E., Guijarro, B., Massutí, E., Alemany, F., Jadaud, A., Perez, J.L., Reglero, P., 2019.  Accounting for ocean connectivity and hydroclimate in fish recruitment fluctuations within transboundary metapopulations.  Ecological Applications 29, e01913.
  * Hilborn (2017) Hilborn, R., 2017.  Are MPAs effective?  ICES Journal of Marine Science 75, 1160–1162.
  * Hutton and Leader-Williams (2003) Hutton, J.M., Leader-Williams, N., 2003\.  Sustainable use and incentive-driven conservation: realigning human and conservation interests.  Oryx 37, 215–226.
  * van Jaarsveld et al. (2005) van Jaarsveld, A.S., Biggs, R., Scholes, R.J., Bohensky, E., Reyers, B., Lynam, T., Musvoto, C., Fabricius, C., 2005\.  Measuring conditions and trends in ecosystem services at multiple scales: the Southern African Millennium Ecosystem Assessment (SAfMA) experience.  Philosophical Transactions of the Royal Society B: Biological Sciences 360, 425–441.
  * Jantke et al. (2018) Jantke, K., Jones, K.R., Allan, J.R., Chauvenet, A.L.M., Watson, J.E.M., Possingham, H.P., 2018\.  Poor ecological representation by an expensive reserve system: Evaluating 35 years of marine protected area expansion.  Conservation Letters 11, e12584.
  * de Juan et al. (2017) de Juan, S., Gelcich, S., Fernández, M., 2017.  Integrating stakeholder perceptions and preferences on ecosystem services in the management of coastal areas.  Ocean & Coastal Management 136, 38–48.
  * de Juan et al. (2015) de Juan, S., Gelcich, S., Ospina-Alvarez, A., Pérez-Matus, A., Fernández, M., 2015.  Applying an ecosystem service approach to unravel links between ecosystems and society in the coast of central Chile.  Science of the Total Environment 533, 122–132.
  * Krueck et al. (2017) Krueck, N.C., Ahmadia, G.N., Green, A., Jones, G.P., Possingham, H.P., Riginos, C., Treml, E.A., Mumby, P.J., 2017\.  Incorporating larval dispersal into MPA design for both conservation and fisheries.  Ecological Applications 27, 925–941.
  * Leslie (2005) Leslie, H.M., 2005.  A Synthesis of Marine Conservation Planning Approaches.  Conservation Biology 19, 1701–1713.
  * Lett et al. (2008) Lett, C., Verley, P., Mullon, C., Parada, C., Brochier, T., Penven, P., Blanke, B., 2008.  A Lagrangian tool for modelling ichthyoplankton dynamics.  Environmental Modelling and Software 23, 1210–1214.
  * Magris et al. (2015) Magris, R.A., Treml, E.A., Pressey, R.L., Weeks, R., 2015\.  Integrating multiple species connectivity and habitat quality into conservation planning for coral reefs.  Ecography 38, 001–016.
  * Martín-López et al. (2012) Martín-López, B., Iniesta-Arandia, I., García-Llorente, M., Palomo, I., Casado-Arzuaga, I., Del Amo, D.G., Gómez-Baggethun, E., Oteros-Rozas, E., Palacios-Agundez, I., Willaarts, B., González, J.A., Santos-Martín, F., Onaindia, M., López-Santiago, C., Montes, C., 2012.  Uncovering Ecosystem Service Bundles through Social Preferences.  PLoS One 7, e38970.
  * Maxwell et al. (2020) Maxwell, S.M., Gjerde, K.M., Conners, M.G., Crowder, L.B., 2020\.  Mobile protected areas for biodiversity on the high seas.  Science 367, 252–254.
  * Millennium Ecosystem Assessment (2005) Millennium Ecosystem Assessment, 2005.  Ecosystems and Human Well-being: Synthesis.  Island Press, Washington, DC.
  * Ministerio del Medio Ambiente Chile (2018) Ministerio del Medio Ambiente Chile, 2018.  Decreto 14 - Aprueba estrategia nacional de biodiversidad 2017 - 2030.  Biblioteca del Congreso Nacional de Chile , 2\.
  * Mora et al. (2006) Mora, C., Andrèfouët, S., Costello, M.J., Kranenburg, C., Rollo, A., Veron, J., Gaston, K.J., Myers, R.A., 2006\.  Coral Reefs and the Global Network of Marine Protected Areas.  Science 312, 1750–1751.
  * Norström et al. (2020) Norström, A.V., Cvitanovic, C., Löf, M.F., West, S., Wyborn, C., Balvanera, P., Bednarek, A.T., Bennett, E.M., Biggs, R., de Bremond, A., Campbell, B.M., Canadell, J.G., Carpenter, S.R., Folke, C., Fulton, E.A., Gaffney, O., Gelcich, S., Jouffray, J.B., Leach, M., Le Tissier, M., Martín-López, B., Louder, E., Loutre, M.F., Meadow, A.M., Nagendra, H., Payne, D., Peterson, G.D., Reyers, B., Scholes, R., Speranza, C.I., Spierenburg, M., Stafford-Smith, M., Tengö, M., van der Hel, S., van Putten, I., Österblom, H., 2020\.  Principles for knowledge co-production in sustainability research.  Nature Sustainability 461, 1–9.
  * O’Neill et al. (2014) O’Neill, B.C., Kriegler, E., Riahi, K., Ebi, K.L., Hallegatte, S., Carter, T.R., Mathur, R., van Vuuren, D.P., 2014\.  A new scenario framework for climate change research: the concept of shared socioeconomic pathways.  Climatic change 122, 387–400.
  * Ospina-Alvarez et al. (2015) Ospina-Alvarez, A., Catalán, I.A., Bernal, M., Roos, D., Palomera, I., 2015.  From egg production to recruits: Connectivity and inter-annual variability in the recruitment patterns of European anchovy in the northwestern Mediterranean.  Progress in Oceanography 138, 431–447.
  * Ospina-Alvarez et al. (2018) Ospina-Alvarez, A., Weidberg, N., Aiken, C.M., Navarrete, S.A., 2018\.  Larval transport in the upwelling ecosystem of central Chile: The effects of vertical migration, developmental time and coastal topography on recruitment.  Progress in Oceanography 168, 82–99.
  * Pérez-Matus et al. (2017) Pérez-Matus, A., Ospina-Alvarez, A., Camus, P.A., Carrasco, S.A., Fernández, M., Gelcich, S., Godoy, N., Patricio Ojeda, F., Pardo, L.M., Rozbaczylo, N., Subida, M.D., Thiel, M., Wieters, E.A., Navarrete, S.A., 2017\.  Temperate rocky subtidal reef community reveals human impacts across the entire food web.  Marine Ecology Progress Series 567, 1–16.
  * Petit et al. (2018) Petit, I.J., Campoy, A.N., Hevia, M.J., Gaymer, C.F., Squeo, F.A., 2018.  Protected areas in Chile: are we managing them?  Revista chilena de historia natural 91, 1.
  * Popova et al. (2019) Popova, E., Vousden, D., Sauer, W.H.H., Mohammed, E.Y., Allain, V., Downey-Breedt, N., Fletcher, R., Gjerde, K.M., Halpin, P.N., Kelly, S., Obura, D., Pecl, G., Roberts, M., Raitsos, D.E., Rogers, A., Samoilys, M., Sumaila, U.R., Tracey, S., Yool, A., 2019.  Ecological connectivity between the areas beyond national jurisdiction and coastal waters: Safeguarding interests of coastal communities in developing countries.  Marine Policy 104, 90–102.
  * Potts et al. (2014) Potts, T., Burdon, D., Jackson, E., Atkins, J., Saunders, J., Hastings, E., Langmead, O., 2014.  Do marine protected areas deliver flows of ecosystem services to support human welfare?  Marine Policy 44, 139–148.
  * Pressey et al. (2007) Pressey, R.L., Cabeza, M., Watts, M.E., Cowling, R.M., Wilson, K.A., 2007.  Conservation planning in a changing world.  Trends in Ecology and Evolution 22, 583–592.
  * Ramesh et al. (2019) Ramesh, N., Rising, J.A., Oremus, K.L., 2019.  The small world of global marine fisheries: The cross-boundary consequences of larval dispersal.  Science 364, 1192–1196.
  * Roberts et al. (2017) Roberts, C.M., O’Leary, B.C., McCauley, D.J., Cury, P.M., Duarte, C.M., Lubchenco, J., Pauly, D., Saenz-Arroyo, A., Sumaila, U.R., Wilson, R.W., Worm, B., Castilla, J.C., 2017\.  Marine reserves can mitigate and promote adaptation to climate change.  Proceedings of the National Academy of Sciences of the United States of America 114, 6167–6175.
  * Rodríguez et al. (2006) Rodríguez, J.P., Beard, T.D., Bennett, E.M., Cumming, G.S., Cork, S.J., Agard, J., Dobson, A.P., Peterson, G.D., 2006\.  Trade-offs across Space, Time, and Ecosystem Services.  Ecology and Society 11.
  * Rossi et al. (2014) Rossi, V., Giacomi, E.S., López, C., García, E.H., 2014\.  Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves.  Geophysical research letters 41, 2883–2891.
  * Saunders et al. (2016) Saunders, M.I., Brown, C.J., Foley, M.M., Febria, C.M., Albright, R., Mehling, M.G., Kavanaugh, M.T., Burfeind, D.D., 2016\.  Human impacts on connectivity in marine and freshwater ecosystems assessed using graph theory: a review.  Marine and Freshwater Research 67, 277–290.
  * Smith and Metaxas (2018) Smith, J., Metaxas, A., 2018\.  A decision tree that can address connectivity in the design of Marine Protected Area Networks (MPAn).  Marine Policy 88, 269–278.
  * Strub et al. (1998) Strub, P.T., Mesías, J.M., Montecino, V., Rutllant, J., Salinas, S., 1998.  Coastal ocean circulation off western South America, in: Robinson, A.R., Brink, K.H. (Eds.), The Sea, 11. John Wiley & Sons, New York, pp. 273–313.
  * Suman et al. (1999) Suman, D., Shivlani, M., Milon, J.W., 1999.  Perceptions and attitudes regarding marine reserves: a comparison of stakeholder groups in the Florida Keys National Marine Sanctuary.  Ocean & Coastal Management 42, 1019–1040.
  * Thiel et al. (2007) Thiel, M., Macaya, E.C., Acuna, E., Arntz, W.E., Bastias, H., Brokordt, K., Camus, P.A., Castilla, J.C., Castro, L.R., Cortes, M., Dumont, C.P., Escribano, R., Fernández, M., Gajardo, J.A., Gaymer, C.F., Gomez, I., Gonzalez, A.E., Gonzalez, H.E., Haye, P.A., Illanes, J.E., Iriarte, J.L., Lancellotti, D.A., Luna-Jorquerai, G., Luxoroi, C., ManrÍquez, P.H., Marin, V., Munoz, P., Navarrete, S.A., Perez, E., Poulin, E., Sellanes, J., Sepulveda, H.H., Stotz, W., Tala, F., Thomas, A., Vargas, C.A., Vasquez, J.A., Vega, J.M.A., 2007\.  The Humboldt Current System of northern and central Chile.  Oceanography and marine biology: an annual review, Vol 35 45, 195–344.
  * UNEP (2010) UNEP, 2010.  Strategic plan for biodiversity 2011-2020 .  Convention on biological diversity , 31\.
  * Visconti et al. (2019) Visconti, P., Butchart, S.H.M., Brooks, T.M., Langhammer, P.F., Marnewick, D., Vergara, S., Yanosky, A., Watson, J.E.M., 2019\.  Protected area targets post-2020.  Science 364, 239–241.
  * Watson et al. (2012) Watson, J.R., Kendall, B.E., Siegel, D.A., Mitarai, S., 2012\.  Changing seascapes, stochastic connectivity, and marine metapopulation dynamics.  The American Naturalist 180, 99–112.
  * Watson et al. (2010) Watson, J.R., Mitarai, S., Siegel, D.A., Caselle, J.E., Dong, C., McWilliams, J.C., 2010\.  Realized and potential larval connectivity in the Southern California Bight.  Marine Ecology Progress Series 401, 31–48.
  * Watts et al. (2009) Watts, M.E., Ball, I.R., Stewart, R.S., Klein, C.J., Wilson, K., Steinback, C., Lourival, R., Kircher, L., Possingham, H.P., 2009.  Marxan with Zones: Software for optimal conservation based land- and sea-use zoning.  Environmental Modelling and Software 24, 1513–1521.
  * White et al. (2014) White, J.W., Schroeger, J., Drake, P.T., Edwards, C.A., 2014\.  The Value of Larval Connectivity Information in the Static Optimization of Marine Reserve Design.  Conservation Letters 7, 533–544.
  * Worm et al. (2006) Worm, B., Barbier, E.B., Beaumont, N., Duffy, J.E., Folke, C., Halpern, B.S., Jackson, J.B.C., Lotze, H.K., Micheli, F., Palumbi, S.R., Sala, E., Selkoe, K.A., Stachowicz, J.J., Watson, R., 2006\.  Impacts of Biodiversity Loss on Ocean Ecosystem Services.  Science 314, 787–790.


# On the approximation of the Hardy $Z$-function via high-order sections

Yochay Jerby  Yochay Jerby, Faculty of Sciences, Holon Institute of
Technology, Holon, 5810201, Israel<EMAIL_ADDRESS>
###### Abstract.

Sections of the Hardy $Z$-function are given by
$Z_{N}(t):=\sum_{k=1}^{N}\frac{cos(\theta(t)-ln(k)t)}{\sqrt{k}}$ for any
$N\in\mathbb{N}$. Sections approximate the Hardy $Z$-function in two ways: (a)
$2Z_{\widetilde{N}(t)}(t)$ is the Hardy-Littlewood approximate functional
equation (AFE) approximation for
$\widetilde{N}(t)=\left[\sqrt{\frac{t}{2\pi}}\right]$. (b) $Z_{N(t)}(t)$ is
Spira’s approximation for $N(t)=\left[\frac{t}{2}\right]$. Spira conjectured,
based on experimental observations, that, contrary to the classical
approximation $(a)$, approximation (b) satisfies the Riemann Hypothesis (RH)
in the sense that all of its zeros are real. We present theoretical
justification for Spira’s conjecture, via new techniques of acceleration of
series, showing that it is essentially equivalent to RH itself.

##  1\. Introduction

###  1.1. Riemann’s Analytic Extension of Zeta

The Riemann zeta function is given by $\zeta(s):=\sum_{n=1}^{\infty}n^{-s}$ in
the range $Re(s)>1$. In his revolutionary 1859 work, Riemann extended the zeta
function $\zeta(s)$ analytically to a meromorphic function on the entire
complex plane with a single pole at $s=1$. This extension allowed him to
explore the zeta function from a complex analytic perspective, uncovering its
profound connection to the distribution of prime numbers. Riemann’s analytic
continuation of $\zeta(s)$ is given by the integral representation:

$\zeta(s)=\frac{\Gamma(1-s)}{2\pi
i}\int_{-\infty}^{\infty}\frac{(-x)^{s}}{e^{x}-1}\frac{dx}{x}.$ (1)

Although the integral representation outlined in equation (1) was of
tremendous importance for Riemman’s theoretical explorations in his
manuscript, it is not amenable for direct computations. To practically compute
$\zeta(s)$, asymptotic techniques are required.

###  1.2. The Approximate Functional Equation (AFE) and the Riemann-Siegel
Formula

The development of the AFE for the $Z$-function due to Hardy and Littlewood,
marks a cornerstone in computation methods of zeta [4, 5, 6, 7]. The Hardy
$Z$-function, denoted as $Z(t)$, is the real function defined by

$Z(t)=e^{i\theta(t)}\zeta\left(\frac{1}{2}+it\right)$ (2)

where $\theta(t)$ is the Riemann-Siegel $\theta$-function, given by the
equation

$\theta(t)=\text{arg}\left(\Gamma\left(\frac{1}{4}+\frac{it}{2}\right)\right)-\frac{t}{2}\log(t),$
(3)

see [2, 8]. The Hardy-Littlewood approximation is encapsulated in the formula

$Z(t)=2\sum_{k=1}^{\widetilde{N}(t)}\frac{cos(\theta(t)-ln(k)t)}{\sqrt{k}}+R(t),$
(4)

where $\widetilde{N}(t)=\left[\sqrt{\frac{t}{2\pi}}\right]$ and the error term
is given by

$R(t)=O\left(\frac{1}{\sqrt{t}}\right).$ (5)

This representation, while powerful, for many purposes, requires further
refinement of the error term $R(t)$ to enhance the level of precision.

In the 1930’s, C. L. Siegel uncovered previously unpublished notes by Riemann
which revealed that, remarkably, Riemann not only knew the Hardy-Littlewood
formula (4) but also developed complex saddle point techniques which enable
the required further evaluation of its error term $R(t)$ for any order [2,
17]. For instance, expanded to first-order, the Riemann-Siegel formula gives

$R(t)=(-1)^{N(t)-1}\left(\frac{t}{2\pi}\right)^{-\frac{1}{4}}\frac{cos(2\pi(p^{2}-p-\frac{1}{16}))}{cos(2\pi
p)}+O\left(\frac{1}{t^{\frac{3}{4}}}\right),$ (6)

where $p=\sqrt{\frac{t}{2\pi}}-\widetilde{N}(t)$. As an asymptotic formula,
for a given specific value of $t$, indefinitely increasing the order of the
expansion does not guarantee improved approximation. It is conjectured by
Berry and Keating [1] that, when expanded to its optimal order, the Riemann-
Siegel formula can reach accuracy level of exponentially decaying error
$O\left(e^{-\pi t}\right)$.

Since its introduction, the Riemann-Siegel formula, especially when enhanced
with the Odlyzko-Schönhage algorithm [13], has been the main method for
numerical verification of the RH, see [3, 11, 14, 12, 16, 20, 21]. The most
comprehensive verification to date, achieved by Platt and Trudigian, confirms
the RH up to $3\cdot 10^{12}$, see [15]. However, it can be argued that one of
the great challenges in the Riemann hypothesis is the fact that while the
Riemann-Siegel formula gives efficient methods for the numerical estimation of
$R(t)$ to any given order, it offers little analytical insight on its
structure, crucial for the understanding of the properties of the zeros in
general, especially as the order of its expansion increases due to the
complexity of the expressions involved.

###  1.3. Sections of the $Z$-function and Spira’s Approximation

In [19, 18] Spira introduced the notion of sections of the AFE of $Z(t)$ which
are defined by

$Z_{N}(t):=\sum_{k=1}^{N}\frac{cos(\theta(t)-ln(k)t)}{\sqrt{k}},$ (7)

and conducted a study of their zeros for any $N\in\mathbb{N}$. In particular,
Spira noted that sections of $Z(t)$ give an additional approximation of the
function to that of the AFE

$Z(t)=Z_{N(t)}(t)+O\left(\frac{1}{\sqrt{t}}\right),$ (8)

at the higher range $N(t)=\left[\frac{t}{2}\right]$. This approximation and
its proof most certainly pre-dates Spira and was probably already known to
Riemann, see Theorem 1.8 of [9] for a variant. A particularly noteworthy
aspect of this additional approximation (8), however, is its unique properties
when compared to the classical AFE, discovered by Spira in his experimental
investigations. For instance, Fig. 1 shows the values of the sections
$Z_{N}(t)$ for the fixed value of $t=3000$ and $N=1,...,1500$ (blue), the
Hardy-Littlewood approximation of $\frac{1}{2}Z(3000)$ (green) and Spira
approximation $Z(3000)$ (brown):

Figure 1. Values of the sections $Z_{N}(t)$ for the fixed value of $t=3000$
and $N=1,...,1500$ (blue), the Hardy-Littlewood approximation of
$\frac{1}{2}Z(3000)$ (green) and the Spira approximation $Z(3000)$ (brown).

For numerical computations, employing the Hardy-Littlewood main sum coupled
with the Riemann-Siegel formula for evaluating $R(t)$ clearly seems
significantly more efficient than utilizing Spira’s approximation. For
instance, for the following reasons:

  1. (1)

The number of terms required for Spira’s approximation,
$N(t)=\left[\frac{t}{2}\right]$, increases quadratically compared to
$\widetilde{N}(t)=\left[\sqrt{\frac{t}{2\pi}}\right]$ required by the Hardy-
Littlewood method, making it far more costly for numerical calculations. This
characteristic alone likely made the approximation (8), even if folklorically
known, impractical for computational use, particularly before the advent of
computers.

  2. (2)

The theoretical asymptotic estimation of Spira’s approximation error, obtained
by classical means, is $O\left(\frac{1}{\sqrt{t}}\right)$, which similar to
that of the Hardy-Littlewood formula, before the application of the Riemann-
Siegel expansion of the error whose first term already gives an error of order
$O\left(t^{-\frac{3}{2}}\right)$. It should be noted, however, that these are
approximate results, and hence this seemingly superior asymptotic bound does
not necessarily ensure greater practical accuracy over Spira’s approximation.

###  1.4. Spira’s conjecture and the Absence of Theoretical Justification

Despite being far more costly, Spira, through his empirical investigations,
discovered that approximation (8) might have a critical significance, see S.6
of [19]. Remarkably, contrary to the Hardy-Littlewood formula, this higher-
range approximation does not seem to admit zeros off the real line. That is,
it is sensitive enough to observe the RH without the need for the expansion of
the error term, as in the classical Riemann-Siegel formula. This observation,
although not explicitly stated by Spira in the following form, can be
formalized as follows:

######  Conjecture (Spira’s RH for sections).

All the non-trivial zeros of $Z_{N(t)}(t)$ are real.

Figure 2 shows $ln\lvert Z(t)\rvert$ (orange) and comparing between the Spira
approximation $ln\lvert Z_{N}(t)\rvert$ with $N=N(t)=205$ (blue - left) and
the classical Hardy-Littlewood approximation $ln\lvert 2Z_{N}(t)\rvert$ with
$N=\widetilde{N}(t)=8$ (blue - right) in the range $412<t<419$:

Figure 2. Graphs of $ln\lvert Z(t)\rvert$ (orange) and the Spira approximation
$ln\lvert Z_{N}(t)\rvert$ with $N=N(t)=205$ (blue - left) and the classical
Hardy-Littlewood approximation $ln\lvert 2Z_{N}(t)\rvert$ with
$N=\widetilde{N}(t)=8$ (blue - right) in the range $412<t<419$

Figure 2 illustrates an instance of two real zeros accurately predicted by
Spira’s formula, as anticipated, but missed by the Hardy-Littlewood formula
(highlighted here in red), which actually erroneously identifies them as
complex zeros

Spira’s conjecture, while grounded in empirical and experimental observations,
lacks direct theoretical underpinnings. Our aim in this work is to present a
theoretical justification for the phenomena observed in Spira’s conjecture of
the advantage of the higher order section $Z_{N(t)}(t)$ over the Hardy-
Littlewood sections $Z_{\widetilde{N}(t)}(t)$, via our new techniques of
accelerated approximations.

##  2\. Spira’s Approximation and Accelerated Approximations

In [10] we have developed an AFE based on the accelerated global series for
$\zeta(s)$ due to Hasse-Sondow, which admits an error of exponential decay.
For the $Z$-function we have the following variant of our accelerated formula:

$Z(t)=\widetilde{Z}_{N(t)}(t)+O\left(e^{-\omega t}\right),$ (9)

where

$\widetilde{Z}_{N}(t):=\sum_{n=0}^{N}\frac{1}{2^{n+1}}\sum_{k=0}^{n}\binom{n}{k}\frac{cos(\theta(t)-ln(k+1)t)}{\sqrt{k+1}},$
(10)

with $\omega>0$ a certain positive constant computed in [10] and
$N(t)=\left[\frac{t}{2}\right]$ is the same order as in Spira’s approximation.
The following Fig. 3 shows $ln\lvert Z(t)\rvert$ (orange) and our
approximation $ln\lvert\widetilde{Z}_{N(t)}(t)\rvert$ with $N=N(t)=205$ (blue)
in the range $412<t<419$.:

Figure 3. Graphs of $ln\lvert Z(t)\rvert$ (orange) and our accelerated
approximation $ln\lvert\widetilde{Z}_{N(t)}(t)\rvert$ of [10] with
$N=N(t)=205$ (blue) in the range $412<t<419$.

Note that our approximation via the accelerated $\widetilde{Z}_{N(t)}(t)$
achieves the superior accuracy expected to be attained by the Hardy-Littlewood
AFE coupled with the Riemann-Siegel formula when evaluated at the optimal
order, according to the Berry and Keating conjecture. Our main result is the
following reformulation:

######  Proposition 2.1.

The following holds

$Z(t)=\widetilde{Z}_{N(t)}(t)+O\left(e^{-\omega t}\right),$ (11)

where

$\widetilde{Z}_{N(t)}(t)=\sum_{k=1}^{N(t)}\widetilde{\alpha}^{acc}_{k}(t)\frac{\cos(\theta(t)-\ln(k)t)}{\sqrt{k}}$
(12)

is the accelerated $N$-th section with the coefficients $\alpha^{acc}_{k}(t)$
given by

$\widetilde{\alpha}^{acc}_{k}(t):=\sum_{n=k-1}^{N(t)}\frac{1}{2^{n+1}}\binom{n}{k-1}.$
(13)

###### Proof.

Denote by

$\beta_{n,k}(t):=\beta^{0}_{n,k}\frac{cos(\theta(t)-ln(k+1)t)}{\sqrt{k+1}},$
(14)

where

$\beta^{0}_{n,k}:=\frac{1}{2^{n+1}}\binom{n}{k}$ (15)

The accelerated formula is given as the sum of all the coefficients
$\beta_{n,k}(t)$ within the triangle $0\leq n\leq N(t)$ and $0\leq k\leq n$.
In particular, summing first along the $k$ indices leads to the definition of

$A(n,t)=\sum_{k=0}^{n}\beta_{n,k}(t).$ (16)

The accelerated formula (9) is thus given by the summation of $A(n,t)$ for
$0\leq n\leq N(t)$. Figure 4 illustrates the triangle of coefficients
$\beta_{n,k}(t)$ and their possible summation orders:

Figure 4. The triangle of coefficients $\beta_{n,k}(t)$ with $A(n,t)$ arising
from the horizontal order of summation (red) and $\widetilde{\alpha}_{k}(t)$
arising from the vertical order of summation (blue).

On the other hand, changing the order of summation to be first along $k-1\leq
n\leq N(t)$ for given $k$, and then along $1\leq k\leq N(t)$ leads to the
definition of $\widetilde{\alpha}_{k}(t)$ and the required formula. ∎

Note that Spira’s section can be written as

$Z_{N(t)}(t)=\sum_{k=1}^{N(t)}\alpha^{step}_{k}(t)\frac{\cos(\theta(t)-\ln(k)t)}{\sqrt{k}},$
(17)

where the coefficients $\alpha^{step}_{k}(t)$ are given by the step function:

$\alpha^{step}_{k}(t):=\left\\{\begin{array}[]{cc}1&1\leq k\leq N(t)\\\
0&\textrm{otherwise}\end{array}\right..$ (18)

Figure 5 shows a comparison between the accelerated coefficients
$\widetilde{\alpha}^{acc}_{k}(t)$ (blue) and the step coefficients
$\alpha^{step}_{k}(t)$ (orange) of Spira’s section for $t=400$ and
$k=1,...,400$:

Figure 5. The accelerated coefficients $\widetilde{\alpha}^{acc}_{k}(t)$
(blue) and the step coefficients $\alpha^{step}_{k}(t)$ (orange) of Spira’s
section for $t=400$ and $k=1,...,400$.

We thus have:

######  Corollary 2.2.

The accelerated section $\widetilde{Z}_{N(t)}(t)$ is related to Spira’s
section $Z_{N(t)}(t)$ via a smoothing of the step coefficients
$\alpha^{step}_{k}(t)$ to obtain $\widetilde{\alpha}^{acc}_{k}(t)$. Moreover,

$lim_{t\rightarrow\infty}\lvert\lvert\widetilde{\alpha}^{acc}_{k}(t)-\alpha^{step}_{k}(t)\rvert\rvert=0$
(19)

in the $\ell_{2}$ norm.

##  3\. Discussion and Conclusions

Sections of the Hardy $Z$-function are known to approximate the Hardy
$Z$-function in two different ways: (a) The classical Hardy-Littlewood AFE
given by $2Z_{\widetilde{N}(t)}(t)$ in (4) and (b) Spira’s approximation given
by $Z_{N(t)}(t)$ in (8). Although the AFE is more efficient and theoretically
expected to be superior, Spira observed through his numerical experimentations
that in practice his approximation consistently satisfies the RH for any
verified $t\in\mathbb{R}$, in the sense that all of its zeros are real. This
is in striking contrast to the properties of the AFE, which requires the
further evaluation of the Riemann-Siegel formula for this purpose.

In this work we presented theoretical justification for this mysterious
observation, based on our new techniques of the asymptotic analysis of series
acceleration, developed in [10] where we introduced a new accelerated
approximation satisfying

$Z(t)=\widetilde{Z}_{N(t)}(t)+O(e^{-\omega t}),$ (20)

where $\omega>0$ is a certain positive constant. This approximation achieves
the superior accuracy expected to be attained by the Hardy-Littlewood AFE
coupled with the Riemann-Siegel formula when evaluated at the optimal order,
according to the Berry and Keating conjecture. In general, for any sequence
$\alpha_{k}=(\alpha_{1},\alpha_{2},...)\in\ell_{2}$ set

$Z(t;\alpha_{k}):=\sum_{k=1}^{\infty}\alpha_{k}\frac{cos(\theta(t)-ln(k)t)}{\sqrt{k}}.$
(21)

In a unified manner, we have shown in Proposition 2.1 that both Spira’s
sections and our accelerated sections can be expressed as elements in this
space

$\begin{array}[]{ccc}\widetilde{Z}_{N(t)}(t)=Z(t;\widetilde{\alpha}_{k}^{acc}(t))&;Z_{N(t)}(t)=Z(t;\alpha_{k}^{step}(t)),\end{array}$
(22)

where $\alpha_{k}^{step}(t)$ is defined in (18) and the definition of
$\alpha_{k}^{acc}(t)$ is defined in (13). Furthermore, Corollary 2.2 shows
that

$lim_{t\rightarrow\infty}\lvert\lvert\widetilde{\alpha}^{acc}_{k}(t)-\alpha^{step}_{k}(t)\rvert\rvert=0.$
(23)

This essentially implies that Spira’s section $Z_{N(t)}(t)$ asymptotically
coincides with our accelerated sections $\widetilde{Z}_{N(t)}(t)$ as
$t\rightarrow\infty$.

In conclusion, our results provide the required theoretical justification to
the numerical phenomena observed by Spira. Indeed, the fact that the
difference between the accelerated sections $\widetilde{Z}_{N(t)}$ and the
Hardy $Z$-function $Z(t)$ itself is imperceptible, coupled with the
convergence of $Z_{N(t)}(t)$ to $\widetilde{Z}_{N(t)}$ at large $t$,
establishes that Spira’s RH for sections is essentially equivalent to RH for
$Z(t)$ itself. Moreover, our results further motivate the definition of the
new parametrized space of sections $Z(t;\alpha_{k})$ and the study of the
properties of their zeros with respect to variation of the parameters
$\alpha_{k}\in\ell_{2}$, suggesting a promising direction for further
exploration.

## Declarations

### Funding

No funding was received to assist with the preparation of this manuscript.

### Conflicts of Interest/Competing Interests

The authors declare that they have no conflict of interest or competing
interests relevant to the content of this article.

### Data Availability

The authors declare that the data supporting the findings of this study are
available within the paper or from the corresponding author upon reasonable
request.

## References

  * [1] Berry, M. V., & Keating, J. P. (1999). The Riemann Zeros and Eigenvalue Asymptotics. SIAM Review, 41(2), 236–266. Society for Industrial and Applied Mathematics.
  * [2] Edwards, H. M. (1974). Riemann’s Zeta Function. Academic Press.
  * [3] Gourdon, X. (2004). The $10^{13}$ first zeros of the Riemann Zeta function, and zeros computation at very large height. Available online.
  * [4] Hardy, G. H., & Littlewood, J. E. (1921). The zeros of Riemann’s zeta function on the critical line. Math. Z., 10, 283–317.
  * [5] Hardy, G. H., & Littlewood, J. E. (1923). The Approximate Functional Equation in the Theory of the Zeta Function, with an Application to the Divisor-Problems of Dirichlet and Piltz. Proceedings of the London Mathematical Society, s2-21(1), 39–74.
  * [6] Hardy, G. H., & Littlewood, J. E. (1929). The approximate functional equations for $\zeta(s)$ and $\zeta^{2}(s)$. Proceedings of the London Mathematical Society, Volume s2-29(1), 81–97.
  * [7] Hardy, G. H., & Littlewood, J. E. (1918). Contributions to the Theory of the Riemann Zeta-Function and the Theory of the Distribution of Primes. Acta Mathematica, 41, 119–196.
  * [8] Ivić, A. (2012). The Theory of Hardy’s Z-Function. Cambridge: Cambridge University Press.
  * [9] Ivić, A. (1985). The Riemann Zeta-Function: Theory and Applications. Dover Publications.
  * [10] Jerby, Y. (2021). An approximate functional equation for the Riemann zeta function with exponentially decaying error. Journal of Approximation Theory, 265, https://doi.org/10.1016/j.jat.2021.105551.
  * [11] Lehman, R. S. (1970). On the Distribution of Zeros of the Riemann Zeta Function. Proc. London Math. Soc., 20, 303–320.
  * [12] Odlyzko, A. M. (1992). The $10^{20}$-th zero of the Riemann zeta function and 175 million of its neighbors. Unpublished book.
  * [13] Odlyzko, A. M., & Schönhage, A. (1988). Fast algorithms for multiple evaluations of the Riemann zeta function. Trans. Amer. Math. Soc., 309(2), 797–809. doi:10.2307/2000939.
  * [14] Odlyzko, A. M. (1987). On the distribution of spacings between zeros of the zeta function. Mathematics of Computation, 48(177), 273–308. doi:10.2307/2007890.
  * [15] Platt, D., & Trudgian, T. (2021). The Riemann hypothesis is true up to $3\cdot 10^{12}$. Bulletin of the London Mathematical Society, Wiley. doi:10.1112/blms.12460.
  * [16] Rosser, J. B., Yohe, J. B., & Schoenfeld, L. (1969). Rigorous Computation and the Zeros of the Riemann Zeta-Function. In Cong. Proc. Int. Fed. Information Process. (pp. 70–76). Washington, DC: Spartan.
  * [17] Siegel, C. L. (1932). Über Riemanns Nachlaß zur analytischen Zahlentheorie. Quellen Studien zur Geschichte der Math. Astron. und Phys. Abt. B: Studien, 2, 45–80. Also in Gesammelte Abhandlungen, Vol. 1. Berlin: Springer-Verlag, 1966.
  * [18] Spira, R. (1967). Zeros of approximate functional equations. Mathematics of Computations, 21(97), 41–48.
  * [19] Spira, R. (1966). Zeros of sections of the zeta function I. Mathematics of Computations, 18, 542–550.
  * [20] Titchmarsh, E. C. (1935). The Zeros of the Riemann Zeta-Function. Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 151(873), 234–255. doi:10.1098/rspa.1935.0146.
  * [21] Turing, A. M. (1953). Some calculations of the Riemann zeta-function. Proceedings of the London Mathematical Society, Third Series, 3, 99–117. doi:10.1112/plms/s3-3.1.99.