# Availability, outage, and capacity of spatially correlated, Australasian
free-space optical networks

Marcus Birch${}^{\,}$ Advanced Instrumentation Technology Centre, Research
School of Astronomy and Astrophysics, Australian National University,
Canberra, ACT 2611, Australia  Corresponding author<EMAIL_ADDRESS>0000-0002-6215-5900 James R. Beattie${}^{\,}$ Research School of Astronomy and
Astrophysics, Australian National University, Canberra, ACT 2611, Australia
Department of Astronomy and Astrophysics, University of California, Santa
Cruz, 1156 High Street, Santa Cruz, CA 96054, United States of America
Corresponding author<EMAIL_ADDRESS>0000-0001-9199-7771 Francis
Bennet  Advanced Instrumentation Technology Centre, Research School of
Astronomy and Astrophysics, Australian National University, Canberra, ACT
2611, Australia  Nicholas Rattenbury  Te Pūnaha Ātea - Space Institute, The
University of Auckland, Auckland 1010, New Zealand  Michael Copeland  Advanced
Instrumentation Technology Centre, Research School of Astronomy and
Astrophysics, Australian National University, Canberra, ACT 2611, Australia
Tony Travouillon  Advanced Instrumentation Technology Centre, Research School
of Astronomy and Astrophysics, Australian National University, Canberra, ACT
2611, Australia  Kate Ferguson  ANU Institute for Space, Australian National
University, Canberra, ACT 2601, Australia  John Cater  Aerospace Research
Centre, University of Canterbury, Christchurch 8140, New Zealand  Mikhael
Sayat  Te Pūnaha Ātea - Space Institute, The University of Auckland, Auckland
1010, New Zealand

###### Abstract

Network capacity and reliability for free space optical communication (FSOC)
is strongly driven by ground station availability, dominated by local cloud
cover causing an outage, and how availability relations between stations
produce network diversity. We combine remote sensing data and novel methods to
provide a generalised framework for assessing and optimising optical ground
station networks. This work is guided by an example network of eight
Australian and New Zealand optical communication ground stations which would
span approximately $60^{\circ}$ in longitude and $20^{\circ}$ in latitude.
Utilising time-dependent cloud cover data from five satellites, we present a
detailed analysis determining the availability and diversity of the network,
finding the Australasian region is well-suited for an optical network with a
69% average site availability and low spatial cloud cover correlations.
Employing methods from computational neuroscience, we provide a Monte Carlo
method for sampling the joint probability distribution of site availabilities
for an arbitrarily sized and point-wise correlated network of ground stations.
Furthermore, we develop a general heuristic for site selection under
availability and correlation optimisations, and combine this with orbital
propagation simulations to compare the data capacity between optimised
networks and the example network. We show that the example network may be
capable of providing tens of terabits per day to a LEO satellite, and up to
99.97% reliability to GEO satellites. We therefore use the Australasian region
to demonstrate novel, generalised tools for assessing and optimising FSOC
ground station networks, and additionally, the suitability of the region for
hosting such a network.

††journal: jocn

##  1 Introduction

Figure 1: Australasian optical ground stations currently being commissioned
(green circles; “the base network"), and proposed for this study (red circles;
“the extended network"), comprising the Australasian Optical Ground Station
Network (AOGSN). Throughout our analysis we use both of these networks as case
studies for availability, correlation, and network optimisation in the
Australasian region.

Australia and New Zealand are pursuing a network of optical ground stations
(OGS) to support the future of optical communication in the region [1]. Free
space optical communication (FSOC) offers myriad advantages over existing
radio frequency technologies, but outages due to clouds are a long-
acknowledged barrier to its widespread adoption and reliability [2]. For
satellite operators seeking to transfer large quantities of data in a timely
manner, an extensive and reliable ground segment is critical. The typical
solution to cloud-induced outages is operating a network of OGS and building
redundancy into the system by providing site diversity such that it is
improbable all OGS are obscured by cloud [3, 4, 5, 6]. In this paper we aim to
develop a means of evaluating the capability of an FSOC network, including
availability, diversity, and ultimately data capacity, alongside a general
framework for network optimisation.

A well-optimised FSOC network is one with high availability at each site, and
minimal correlation of cloud cover between sites for maximising diversity. Low
Earth orbit (LEO) satellites can improve data throughput by engaging with a
well-optimised network, increasing the number of possible links. For
Geostationary orbit (GEO) satellites and deep space, where mutual visibility
can be simultaneously maintained to multiple OGS, a well-optimised network
boosts reliability, i.e. the probability that at least one OGS is available.

We approach this issue by using the proposed Australasian Optical Ground
Station Network (AOGSN) as an example [1]. The AOGSN will include at least the
Australian National University’s Optical Communication Ground Station at Mount
Stromlo Observatory, the Australian Defence Science & Technology Group OGS in
Edinburgh, South Australia, and the University of Western Australia’s OGS in
Perth [1, 7, 8, 9]. All three of these are in a state of commissioning and aim
to enter operation in 2023. To create a robust model network to study, we
propose an additional five sites at Mt Kent Observatory in Queensland, Mt John
Observatory in New Zealand, Alice Springs in the Northern territory, Mildura
in Victoria, and Dongara in central Western Australia. These additional sites
are chosen as places with existing observatories or satellite ground segment
infrastructure, and represent realistic OGS placements. The base and extended
AOGSN will be used henceforth for reference to the initial three OGS as a
network or all eight OGS respectively. Figure 1 shows the sites of this
proposed network. We use the example of the extended AOGSN throughout this
paper to compare with a network chosen from an optimisation model, and to test
methods for availability, diversity, and capacity estimation.

Complex theoretical models have been developed for determining and predicting
outages, including factors such as fading of the atmospheric channel due to
turbulence [3]. An experimental campaign to measure atmospheric turbulence at
selected sites is considered in Section 8, however this paper focuses on a
spatial availability model without prior site selection. Therefore, as with
earlier literature, we adopt a simple model for site outages: An OGS
experiences an outage only due to cloud cover, which prevents the ground
station from linking with the orbiting satellite [10, 11].

Unlike clouds, where associated outages are a generic issue for any network,
turbulence-induced outages are strongly system-specific. They may depend upon
the orbital height of the satellite, details of the telescope and ground
station specifications, e.g adaptive optics and uplink or downlink
capabilities, etc [12, 13, 14, 15]. This means that our generic, cloud-based
outage model will require careful comparison with a real system, on a case-by-
case basis. For example, an adaptive optics-compensated downlink from GEO, at
a high quality site, is unlikely to experience many turbulent outages and will
therefore have availability close to the cloud availability. Alternatively,
uplinking to a LEO satellite in the same atmospheric conditions could suffer
turbulent degradation that prevents the AO system from closing the loop for
long periods of time, causing an outage. The difficulties of predicting
turbulence conditions, and how this could be expanded into a rigorous
turbulent and cloud availability model are discussed further in Section 8. Our
current model therefore encapsulates availability as the probability that a
site is not cloudy at any given time. Given atmospheric turbulence might
create additional outages, our model provides an upper bound on a realistic
availability, which may be required for a detailed system-specific network
performance analysis, which is beyond the scope of this study.

We provide a more rigorous definition of this cloud model based on the
Bernoulli distribution in Section 2. Remote sensing satellites measure cloud
cover in different ways and with different issues, prompting us to adopt an
approach that incorporates numerous satellites at different temporal and
spatial resolutions, with decades of data in time. Deriving network diversity
requires knowledge of the point-wise cloud correlation between sites.
Simultaneous cloud cover measurements of the entire Australasian region are
only possible with GEOs given their coverage, therefore this study places
emphasis on data acquired from Himawari-8, a GEO weather satellite operated by
the Japan Meteorological Agency, throughout our analysis [16].

Figure 2: Availability of Australasian region as detected by AHI/Himawari-8,
averaged from using twice daily cloud cover masks from 2015 to 2022. White
dots correspond to geographical position of all OGS of the extended AOGSN, and
green outlines indicate base AOGSN sites. Annotated availabilities are shown
for each site. Data provided by Earth Observation Research Centre (Japan
Aerospace Exploration Agency).

We develop a general method for computing outage probabilities for an
arbitrary sized network with arbitrary point-wise correlations between each
node, analogous to a network of spiking biological neurons with a specified
covariance structure [17]. We sample our model using the availability and
point-wise correlation from the base and extended AOGSN to provide robust
outage statistics, extending the analysis performed in [10]. Furthermore,
utilising our availability and diversity statistics we compare the base and
extended AOGSN outage probabilities with a network that is optimised to
maximise availability and minimise point-wise correlation for every ground
station (node). Diversity optimisation methods have been recently developed
[3, 6, 11, 18]. However, we present a novel model for network optimisation
that is spatially-resolved, not requiring pre-selected site locations.

This study is organised by primarily using the AOGSN model network as an
example to drive methods and provide realistic outputs, and is structured as
follows. In Section 2 we measure cloud cover over the entire Australasian
region, and for sites of the AOGSN. Section 3 covers the use of high temporal
resolution, and simultaneous cloud measurements to directly measure spatial
correlation of cloud cover, and correlations between selected sites. We then
present a novel Monte Carlo method in Section 4 for determining the outage
probabilities of a network based on the availabilities and correlations
between ground stations. A simple heuristic model for site selection to
maximise availability and diversity is presented in Section 5, then used to
estimate an optimised network which is compared with the AOGSN across a number
of metrics. Finally, in Section 6, we use orbital simulations to estimate the
coverage of the network, and statistical link duration per day with LEO
satellites. Estimates of network-wide link duration from simulation are
combined with availabilities to estimate network capacity. For GEOs, we use
our outage probability model to determine the network outage probability as a
function of orbital longitude in Section 7. Finally, in Section 8, we
summarise and list the key results of the study.

##  2 Cloud cover and site availability

Cloud cover, $\mathcal{C}$, is the major factor in degrading an FSOC network’s
availability, $A$ [2, 19]. Although satellites measure $\mathcal{C}$ in
different ways, and many data products we use have already undergone spatial
or time averaging, they are all capable of producing a binary cloud mask which
flags $\mathcal{C}$ for each pixel, $(x_{i},y_{j})$, and time, $t_{n}$, in a
satellite image. We consider that $\mathcal{C}$ is instantaneously measured as
0 (no cloud detection) or 1 (cloud detection) for each $(x_{i},y_{j})$ and
$t_{n}$. Hence, over a time-interval, $t\in[t_{1},t_{n}]$, the cloud cover is
a vector $\bm{C}$ of $\left\\{0,1\right\\}$, with length $n$. Each component
of $\bm{C}=(C_{1},C_{2},\ldots C_{n})$, $\mathcal{C}_{n}$ follows a Bernoulli
distribution, $\mathcal{C}_{n}\sim\text{Bernoulli}(\Omega)$, where
$p(\mathcal{C}_{n}=1)=\Omega$, and $p(\mathcal{C}_{n}=0)=1-\Omega=A$, hence,

$\displaystyle p(\mathcal{C};\,\Omega)=\Omega^{\mathcal{C}}A^{1-\mathcal{C}},$
(1)

is the probability distribution function for $\mathcal{C}$. This is true for
each $(x_{i},y_{j})$ on the grid.

Let us consider a network $\mathcal{N}$ with $N$ nodes or sites. We call the
$k^{\text{th}}$ site $\mathcal{S}_{k}$, and hence
$\mathcal{N}\equiv\left\\{\mathcal{S}_{1},\mathcal{S}_{2},\ldots,\mathcal{S}_{N}\right\\}$
is the set of $N$ sites, each with grid coordinates $(x_{i},y_{j})$. We may
choose to index the sites themselves, rather than the grid coordinates, e.g.,
the probability distribution function at site $k$ is
$p_{k}(\mathcal{C}_{k};\Omega_{k})$111We avoid a subscript of $\mathcal{C}$
for the time vector component index, $n$, which would interfere with our $k$
site index for simplicity.. Because at any point in time each grid coordinate
follows a Bernoulli distribution for $\mathcal{C}$, so does each site.
Therefore, the probability density function for each $k$ site is
$p_{k}(\mathcal{C};\Omega_{k})=\Omega_{k}^{\mathcal{C}_{k}}A_{k}^{1-\mathcal{C}_{k}}$.
Note, directly from Equation 1 we can immediately estimate the availability
$A$ at site $k$ through the time-average222We use the ensemble
$\left\langle{X}\right\rangle_{t}$ notation for the time-average of the
quantity $X$.
$\left\langle{\bm{\mathcal{C}}_{k}}\right\rangle_{t}=\Omega_{k}=1-A_{k}$.

$A$ has the most intuitive interpretation for FSOC ground stations as it
corresponds to the probability a site is available for communications. We
conduct an extensive study of $\mathcal{C}$ from numerous satellites to
determine $A$ for each ground station of the AOGSN. Five separate remote
sensing satellites and instruments are used to determine $A$, and long time
spans are used to minimise climatic variation. Table 1 summarises the
satellites and instruments used for this study and the relevant time span of
data associated with each. We use derived cloud retrieval products, with a
focus on Sun-synchronous LEO satellites, and Himawari-8, which has been the
prime GEO for use over the Australasian region since 2015 [16].

Although Earth observation LEO satellites can produce cloud mask swathes with
$<100$ m spatial resolution, we find working with low resolution global
averages easier given the decades long time spans involved and subsequent
quantity of data. Comparisons with higher resolution data products, e.g. from
Himawari-8, found reasonable agreement, likely from the large time span of
data involved and the relatively homogeneous terrain of Australia (Appendix
B). New Zealand is naturally more impacted by spatial averaging; subsection
2.2 outlines how Aqua/AIRS and MERRA-2 data is discarded for the Mt John node
in the AOGSN. Therefore, for the LEO satellites: Terra+Aqua/MODIS, Aqua/AIRS,
and Suomi-NPP/VIIRS data, we use $1^{\circ}\times 1^{\circ}$ global grids with
monthly $\Omega$ instead of raw $\mathcal{C}$ measurements at high resolution
[20, 21, 22, 23].

Spatial averaging for terrain more varied than the Australian continent is
likely to induce more errors, so doing the full computation on individual
swathes would be preferred. MERRA-2, a modern retrospective analysis of
assimilated satellite data since 1980, is able to provide a uniquely long time
span with daily cloud information at $0.5^{\circ}\times 0.6^{\circ}$ spatial
resolution [24]. However, MERRA-2’s reliance on older satellite data and large
scale spatial averaging and interpolation can produce erroneous results.
MERRA-2 is analysed before including in any statistical means (e.g. for the Mt
Stromlo and Alice Springs station it produces outlier results compared with
other data). Long time span data, such as MERRA-2, is important for reducing
uncertainty from periodic climatological phenomena such as the El Niño
Southern Oscillation and Indian Ocean Dipole [25]. However, anthropogenic
climate change is a non-periodic climatogolical effect that can induce
significant error on long time span data, e.g. MERRA-2, as the effect on
clouds and their circulation has been studied extensively [26, 27, 28, 29].

Table 1: Remote sensing satellites used in this paper, and spatial resolution
of relevant cloud retrieval data products used.

Satellite |  Instrument |  Time span |  Spatial Resolution
---|---|---|---
Himawari-8 |  AHI |  2015- |  5 km
Terra/Aqua |  MODIS |  2002- |  $1^{\circ}$
Aqua |  AIRS |  2003- |  $1^{\circ}$
Suomi NPP |  VIIRS |  2012- |  $1^{\circ}$
MERRA-2 |  Assimilation |  1980- |  $0.5^{\circ}$

  * •

Notes. MERRA-2 is a historical assimilation of data from numerous sources
provided by the NASA Global Modelling and Assimilation Office [24].

###  2.1 Spatially resolved Himawari-8 maps

Himawari-8 is used for establishing an initial, spatially-resolved picture of
$A$ over the Australasian region, given it is a GEO satellite that produces a
cloud mask of the entire region every ten minutes at a 5 km resolution [16].
Due to the large quantity of data produced by Himawari-8 since 2015, images at
9am and 3pm Australian Eastern Standard Time (UTC+10) are used, an Australian
convention established by the Bureau of Meteorology. Figure 2 shows $A$ over
eight years of Himawari-8 daytime cloud product provided by JAXA/EORC333Japan
Aerospace Exploration Agency, Earth Observation Research Centre and the Japan
Meteorological Agency [30, 31]. Cirrus cloud identified by Himawari-8 is
considered to be clear sky for Figure 2 given its low $1.55$ µm attenuation
[32]. Himawari-8 can produce misleading results for pixels with high altitude
areas that typically sit above cloud height, seen for the Alps down the length
of New Zealand’s South Island in Figure 2 [31], but this does not affect any
area of Australia.

###  2.2 Seasonal variation with collated data

Figure 3: Seasonal availability for all ground stations in extended AOGSN. The
availability, $A$, is averaged from satellites tabulated in Table 1. The
$1\sigma$ fluctuations, visualised with the blue band, display the variability
between the different satellites. Uncertainty in the annual average of $A$
combines the $\sigma$ from satellite and the $\sigma$ over the monthly values.
Figure 4: Surfaces of the Pearson correlation coefficient,
$r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})$, between site $k$ in the
network and all other pixels,
$\bm{\mathcal{C}}_{ij}=\bm{\mathcal{C}}(x_{i},y_{j})$, for the Australasian
region. Correlation surfaces are constructed by comparing the time-series of
daily AHI/Himawari-8 cloud retrieval images between 2015-2022. We show
contours for $r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})=0.2$ and
$0.4$, overlaid on each surface.

To provide a more robust estimate of $A$, we then calculate $A$ from all five
satellites and decades of data for all sites of the extended AOGSN. All
sources of remote sensing data outlined in Table 1 were used except for the Mt
John OGS where the Aqua/AIRS and MERRA-2 cloud retrieval algorithms produced
erroneous results, likely due to nearby snow and ice or low resolution pixels
including high altitude mountains. For low spatial resolution images, the
closest pixel to the site is used to represent the sites cloud cover. Given
Himawari-8’s 500 m/pixel spatial resolution, a region-of-interest (ROI) is
drawn around the sites and pixels within that ROI are averaged per frame. We
model the spatial extent of these ROIs to correspond with a projected
horizontal distance from a single layer of cloud at some height when observed
at $\theta_{\text{elevation}}=30^{\circ}$444
$\theta_{\text{elevation}}=30^{\circ}$ is considered the minimum elevation
angle for space-to-ground links throughout this paper.. Himawari-8 is capable
of detecting cloud height, so we measure 5 km as the mean cloud height over
the entire Australasian region from 2015-2022, such that our ROI for
Himawari-8 is approximately 9 km in radius, or $\pm 2$ pixels.

Figure 3 shows $A$ across all satellites and the annual average of
$A\pm\sigma$, with $\sigma$ from satellite and seasonal variation. From Figure
3 we see that most sites have high availability with an average of $A=0.69$,
particularly the arid sites such as Alice Springs ($A=0.83$) and Dongara
($A=0.77$). Most sites exhibit a drop in $A$ for winter months, except for the
northern sites of Alice Springs ($23.8^{\circ}S$) and Mt Kent
($27.8^{\circ}S)$, which experience a partial dry and wet seasonal cycle. Mt
John is an outlier due to its low availability and pattern of clearer skies in
the winter despite its southern latitude ($44^{\circ}$S). Appendix B
demonstrates the independent results of $A$ from the different data sets used.

##  3 Network diversity

###  3.1 Spatial correlation of cloud cover

Correlation of cloud cover between sites is the primary impact on network
diversity, as two or more network nodes may be affected by the same weather.
An effective FSOC network should be designed such that its OGS do not have
strongly correlated cloud cover. For the case of GEO or deep space
communication where visibility is maintained to an entire network, diversity
corresponds directly to the likelihood that a link can be established with
said network.

We use Himwari-8 for all cloud correlation analysis as it has made
simultaneous $\mathcal{C}$ measurements of the entire region over the last
eight years. Each pixel is flagged as $0$ if no cloud or only cirrus cloud is
present, otherwise it is flagged as $1$, as we described in Section 2.
Therefore we have the vector $\bm{\mathcal{C}}_{ij}$555Throughout this study,
we use the subscript $X_{ij}$ for function $X(x_{i},y_{j})$ evaluated at grid
coordinates $(x_{i},y_{j})$. of length $n$, where $n$ is the number of daily
images from 2015-2022, for each $5$ km$\times 5$ km pixel $(x_{i},x_{j})$ in
Himawari-8 frames of Australasia. We calculate the point-wise Pearson (linear)
correlation coefficient, $r_{ij}$, between site $k$, $\bm{\mathcal{C}}_{k}$,
and every other pixel in Australasia, $\bm{\mathcal{C}}_{ij}$.

$r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})=\frac{\left\langle{\Big{(}\bm{\mathcal{C}}_{k}-\Omega_{k}\Big{)}\Big{(}\bm{\mathcal{C}}_{ij}-\Omega_{ij}\Big{)}}\right\rangle_{t}}{\left[\Big{\langle}\Big{(}\bm{\mathcal{C}}_{k}-\Omega_{k}\Big{)}^{2}\Big{\rangle}_{t}\Big{\langle}\Big{(}\bm{\mathcal{C}}_{k}-\Omega_{ij})\Big{)}^{2}\Big{\rangle}_{t}\right]^{1/2}},$
(2)

where $-1\leq r_{ij}\leq 1$. For $r_{ij}=1$ the cloud cover between two pixels
is perfectly linearly correlated666Note that this type of correlation only
captures linear correlations, and more complex, non-linear correlations may
exist but will not be accounted for with this statistic, see e.g. [33]. in
time, hence when one pixel is cloudy, so is the other one. Conversely, when
$r_{ij}=-1$ the pixels have perfect linear negative correlation, and when one
pixel is cloudy the other is not. Finally, $r_{ij}=0$ corresponds to the case
where two pixels are not linearly correlated in time.
$r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})$ is therefore a two-
dimensional surface describing how $\mathcal{C}$ throughout the region is
spatially correlated around site $k$. The length scale and shape of $r_{ij}$
is driven by localised and global weather dynamics, e.g. atmospheric
circulation cells on the global scale and rain shadow more locally.

Figure 4 shows $r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})$ for each
of site in the extended AOGSN. We find contours of equal $r_{ij}$ around sites
that represent the weather dynamics driving $\mathcal{C}$, e.g. elongation
along the East-West axis for southern sites due to prevailing westerlies
across Australia. From Figure 4, we see that Alice Springs experiences
correlated $\mathcal{C}$ over an extremely large area, the result of how
pressure systems and fronts traverse the desert without coastal or topological
effects. It is also evident that Dongara and Perth are located within each
others $r_{ij}=0.4$ contour, implying they will be strongly correlated as we
show in subsection 3.2. The Mt John site highlights the effect of rain shadow,
as the correlation scale is stretched along the length of the New Zealand
Southern Alps. Mt John correlations also show misleading results on the alps,
due to those pixels sitting above cloud height. Negative correlation between
northern Australia and south western Australia is also clear in Figure 4,
which occurs because northern Australia experiences has a tropical
seasonality, i.e. a winter dry season, unlike the rest of the continent.

This method of computing $r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})$
for $k$ sites is an extension earlier work [10], i.e. we compute $r_{ij}$
against $\mathcal{C}_{k}$ for the whole grid $(x_{i},y_{j})$, and not just
point-wise between preselected site locations. This spatially-resolved
approach, despite being more computationally expensive, allows us to construct
a two-dimensional surface of $r_{ij}$ for each site. In Section 5, we show how
this quantity can be used to generate a network of maximally diverse ground
stations.

###  3.2 Site-specific correlations

To probe the pair-wise correlations between sites we calculate
$r(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{l})$ between all sites of the
extended AOGSN, where $k,l\in\\{1,\dots,N\\}$ for a network of $N$ sites. We
show the results with a correlation matrix in Figure 5. Figure 5 shows
typically low, $r\lessapprox 0.1$, between sites. However, we do note the
Mildura-Edinburgh ($r=0.51$) and Perth-Dongara ($r=0.5$) sites show relatively
strong, positive correlations. For example, this corresponds to half of all
outages at Perth or Dongara affecting both sites. Negative correlation between
Alice Springs and Perth/Dongara is due to the winter dry season of northern
Australia, visible in Figure 4.

Figure 5: Pearson correlation matrix (Equation 2 for each pair-wise site) for
all sites in extended AOGSN calculated using daily cloud retrieval frames of
Australasian region since 2015 from Himawari-8.

As a simple metric to understand how diverse a network of $N$ sites is, we
compute the mean absolute correlation in the network,

$\bar{|r|}={N\choose
2}^{-1}\sum_{k=1}^{N-1}\sum_{l=k+1}^{N}\Big{|}r(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{l})\Big{|},$
(3)

introduced earlier [10]777Note, that compared to the value of $\bar{r}$ in
[10], we choose to take the mean of the absolute value, $\bar{|r|}$ which
means we are probing the average magnitude of the correlation over the network
as opposed to the average correlation. The latter may be zero for a network
that has both strong positive and negative correlations, hence $\bar{r}=0$
does not imply that the network is, as a whole, has weak correlations, but
$\bar{|r|}=0$ does.. We find $\bar{|r|}=0.088\pm 0.11$ and $\bar{|r|}=0.073\pm
0.076$ for the extended and base AOGSNs respectively, which is reasonably low
despite being skewed by a number of strongly correlated sites (e.g., Mildura-
Edinburgh) that were shown in Figure 5. Negatively correlated sites contribute
given the absolute approach of Equation 3, and the detrimental effects of
negative correlations between network sites is explored in Appendix C.
Nonetheless, $\bar{|r|}=0.088$ is representative of a well-diversified
network, owing to the large area of the Australasian region. This highlights
the unique position for Australia and New Zealand to offer a robust, highly-
diverse FSOC ground segment, when combined with the high availabilities
estimated in Section 2.

##  4 Network reliability and outage probability

The diversity and availability of a network can be used to estimate the
network reliability and outage probability. To understand the outage
probability, we construct $p_{\text{avail}}(k\leq M)$, the probability that at
least $M$ sites in the network of $N$ nodes ($M\leq N$) are available in a
network (the cumulative distribution function of the network availabilities.).
The case of $M=0$ provides the network-wide outage probability,
$p_{\text{avail}}(M=0)$. For GEOs and deep space targets, which can maintain
visibility with an entire hemisphere $1-p_{\text{avail}}(M=0)$ corresponds
directly to the network reliability. Figure 5 shows that some significant
correlations exist between AOGSN nodes, as do the networks proposed in [10].
Therefore, it is necessary to construct $p_{\text{avail}}(k\leq M)$ including
point-wise correlations between sites.

Constructing a closed-form analytical expression for $p_{\text{avail}}(k\leq
M)$ with correlation maintained between greater than three sites is impossible
[34]. An earlier study provided an expression for uncorrelated sites and
empirically measured the outage probabilities from $\bm{\mathcal{C}}$ data for
the correlated case [10]. Empirically estimating $p_{\text{avail}}(k\leq M)$
can be compromised by the length of time required to determine low probability
events such as a network-wide outages [10]. For our case, we use only twice
daily images, and find only one instance of a network-wide outage over the
2015-2022 time span for the extended AOGSN. This means empirically, similar to
[10], we cannot provide a robust estimate of $p_{\text{avail}}(M=0)$ and
require a better method to estimate these low probability, correlated events.
A robust method can also be used to estimate from correlated weather data
built up over a relatively short amount of time.

Because there is no closed-form analytical expression for the probability
function of a correlated network of Bernoulli events, we opt for a Monte Carlo
sampling approach. Our goal is to construct $p_{\text{avail}}(k\leq M)$ from a
distribution that we are able to randomly draw events from. As per Section 2,
each event is a time realisation of the cloud cover at each site,
$\mathcal{C}_{k}$. We therefore have to construct the the joint availability
distribution (a multivariate Bernoulli distribution of dimension $N$) with
correlations between each $\mathcal{C}_{k}$ marginal distribution. Again,
$k,l\in\\{1,\dots,N\\}$ for a network of $N$ sites. We do this by utilising a
neuronal spike model [17], noting that a ground station network described by a
correlated network of Bernoulli events is directly analogous to the correlated
network of firing neurons [17].

Following [17], the method of constructing and sampling from the joint
availability distribution is as follows. Consider $\bm{x}\in\mathbb{R}^{N}$,
with components $x_{k}$. Let $\bm{x}$ follow a multivariate Gaussian
distribution (that is, each $x_{k},$ admits to a Gaussian distribution),

$\displaystyle
f(\bm{x};\bm{\mu},\bm{\Lambda})=(2\pi)^{-\frac{n}{2}}\text{det}(\bm{\Lambda})^{-\frac{1}{2}}\exp\left\\{-\frac{1}{2}(\mathbf{x}-\bm{\mu})^{T}\bm{\Lambda}^{-1}(\bm{x}-\bm{\mu})\right\\},$
(4)

with mean $\bm{\mu}\in\mathbb{R}^{N}$ and covariance
$\bm{\Lambda}\in\mathbb{R}^{N\times N}$, a symmetric and positive-definite
matrix and $\text{det}(\bm{\Lambda})$ is the matrix determinant. [17] showed
that if we are able to sample $f(\bm{x};\bm{\mu},\bm{\Lambda})$ and pick a
threshold,

$\displaystyle\mathcal{C}_{k}=1,\;\;\text{if}\;\;x_{k}>0,$ (5)

then we are able to simulate a correlated multivariate Bernoulli process, with
$\Omega_{k}$ determined by $\mu_{k}$ and the covariance $\Gamma_{kl}$ between
each ground station determined by $\Lambda_{kl}$. Assuming that
$\Lambda_{kk}=1$, the relations between these quantities are as follows,

$\displaystyle\Omega_{k}$ $\displaystyle=\phi(\mu_{k}),$ (6)
$\displaystyle\Gamma_{kk}$ $\displaystyle=\phi(\mu_{k})\phi(-\mu_{k}),$ (7)
$\displaystyle\Gamma_{kl}$
$\displaystyle=\phi_{2}(\mu_{k},\mu_{l},\Lambda_{kl})-\phi(\mu_{k})\phi(\mu_{l}),\;k\neq
l,$ (8)

where

$\displaystyle\phi(\mu_{k})=\frac{1}{2}\left[1+\text{erf}(\mu_{k}/\sqrt{2})\right],$
(9)

is the cumulative distribution function for a univariate, standardised
($\mu=0$; $\Gamma=1$) Gaussian distribution, and
$\phi_{2}(\mu_{k},\mu_{l},\Sigma_{kl})$ is the two-dimensional counterpart.
$\mu_{k}$ is determined simply by inverting Equation 9, which can then be used
to generate $\Gamma_{kk}$ using Equation 7. For the off-diagonal elements of
$\bm{\Lambda}$, we perform a root-finding algorithm to determine what values
of $\Lambda_{kl}$ satisfy

$\displaystyle\Gamma_{kl}-\left[\phi_{2}(\mu_{k},\mu_{l},\Lambda_{kl})-\phi(\mu_{k})\phi(\mu_{l})\right]=0,$
(10)

noting that $-1\leq\Lambda_{kl}\leq 1$ (this is true for Bernoulli events, see
Appendix A in [17], but not generally true for other random processes) and
that $\Gamma_{kl}$ is determined directly from the data (i.e., we compute the
covariance between ground stations of the network, which from Equation 2 one
can see that it is
$[\langle(\bm{\mathcal{C}}_{k}-\Omega_{k})^{2}\rangle_{t}\langle(\bm{\mathcal{C}}_{l}-\Omega_{l}))^{2}\rangle_{t}]^{1/2}r(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{l})$.
$\Gamma_{kl}$ and $\Lambda_{kl}$ are monotonic functions of one another, and
therefore Equation 10 has a single root between $-1$ and $1$, which we
determine using scipy [35]. This method differs from [5] by numerically
solving Equation 10 instead of obtaining an analytical expression by
approximation.

Figure 6: Cumulative distribution functions $p_{\rm avail}(k\leq M)$, for
$k\leq M$ available sites, constructed using our Monte Carlo method outlined
in Section 4 for the extended (solid black line; $M=8$) and base (solid grey
line; $M=3$) AOGSN, including point-wise correlations between ground stations
(see Section 4). We show uncorrelated ($r=0$; black dashed line) and highly
correlated ($r=0.5$; black dotted line) network availability results to
compare with the extended and base AOGSN. We have highlighted the special case
of $p_{\text{avail}}(M=0)$ (red and green band) representing a network-wide
outage.

Utilising this construction and sampling method, we create three multivariate
availability distributions for the extended AOGSN, and three availability
cumulative distribution functions by sampling the availability distribution
$10^{8}$ times (equivalent to a $10^{8}$ time realisations), all which use
$\Omega_{i}$ directly determined from the data. The difference between the
three distributions is as follows: (1) we construct a distribution with
$\bm{\Gamma}$ directly from our data, encoding the real
$r(\mathcal{C}_{k},\mathcal{C}_{l})$ between sites; (2) we run a
$r(\mathcal{C}_{k},\mathcal{C}_{l})=0,\;k\neq l$ simulation, which enforces
that there is no correlation between sites, comparable to the results in [10];
and (3) we choose a $r(\mathcal{C}_{k},\mathcal{C}_{l})=0.5,\;k\neq l$ to
understand the outage probability for a network that is strongly correlated in
$\mathcal{C}$.

Figure 6 shows the cumulative distribution function $p_{\text{avail}}(k\leq
M)$ for $M=0-8$, with cases (1), (2), and (3) for the extended AOGSN, and with
only case (1) for the base AOGSN. We estimate network-wide outage
probabilities of $p_{\text{avail}}(M=0)=2.4\times 10^{-4}$ and
$p_{\text{avail}}(M=0)=0.064$ for the extended and base AOGSN configurations
respectively. From Figure 6 we see that a small ($N=3$) and well-separated
network, e.g., the base AOGSN, can have (marginally) lower or comparable
outage probability than a highly available but very correlated network of
eight OGS. Figure 6 also shows while $p_{\text{avail}}(M=0)=2.4\times 10^{-4}$
for measured $r(\mathcal{C}_{k},\mathcal{C}_{l})$,
$p_{\text{avail}}(M=0)=6.1\times 10^{-5}$ for case (2), i.e., uncorrelated
sites with $r(\mathcal{C}_{k},\mathcal{C}_{l})=0$. An analytical expression
exists for $r=0$ [10], but we can see that realistically modelling
correlations results in a $4\times$ higher probability of a network-wide
outage. Furthermore, for more correlated and dense networks, e.g. the Japanese
network outlined in [36]. This discrepancy can be orders of magnitude as we
see with case (3) for $r(\mathcal{C}_{k},\mathcal{C}_{l})=0.5$. Appendix C
shows a negatively correlated network, $r(\mathcal{C}_{k},\mathcal{C}_{l})<0$,
demonstrating how any correlations (negative or positive) increase
$p_{\text{avail}}(k\leq M)$ compared to a
$r(\mathcal{C}_{k},\mathcal{C}_{l})=0.0$ network, in agreement with [5].

This method shows how robust estimates can be produced for the network outages
of a correlated network using availability and correlation data. It also
highlights how a sparse, highly available, and diverse network, such as the
example network we study, can provide an extremely reliable ground segment,
overcoming one of the most significant barriers for space-to-ground FSOC.
Next, we consider the case where we want to generate our own optimal network
that maximises availability and diversity.

##  5 Generating an optimal network

We are interested in developing a method for site selection in an FSOC
network. An optimised network, produced without preselection, can also be
compared with an existing network, e.g. our example AOGSN, to judge how well-
optimised the network is. This method could also be adopted to select
additional sites in an established network. Methods have been presented for
how diversity can be optimised with varying channel and outage models, but
their optimisation relies on preselected sites in a network [3, 6].

Figure 4 indicates that correlation surfaces can be used for optimally
planning an FSOC network. We expand this idea by considering a simple
spatially-resolved model for iteratively determining the coordinates
$(x_{i},y_{j})$ of $N$ sites $\mathcal{S}_{k}$, in the network $\mathcal{N}$
(the same notation as we used in Section 2). Our method for determining the
optimal $(x_{i},y_{j})$ for each $\mathcal{S}_{k}$ in the network relies upon
two site characteristics: the first is average cloud cover
$0\leq\Omega_{k}\leq 1$ at each $\mathcal{S}_{k}$. The second is the temporal
correlation $-1\leq r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})\leq 1$.
We discuss the exact procedure for generating the correlation in Section 3. To
maximise the connectivity of each $\mathcal{S}_{k}$ and robustness of
$\mathcal{N}$, we propose a simple heuristic: (1) the optimal site is one that
minimises the amount of cloud coverage, hence maximising the time it may be
available on the $\mathcal{N}$, and (2), the optimal $\mathcal{N}$ is one that
supports a series of uncorrelated $\mathcal{S}_{k}$ (supported by the lowest
probability of an outage in Figure 6), such that if any $\mathcal{S}_{k}$ are
unavailable due to cloud cover events, other sites in $\mathcal{N}$ will most
likely not suffer the same deficiency.

Based on this simple heuristic, the first site $\mathcal{S}_{1}$ we pick is
located at,

$\displaystyle\mathcal{S}_{1}(x_{i},y_{j})=\text{argmin}_{ij}\left\\{\Omega_{ij}\right\\},$
(11)

where $\text{argmin}_{ij}$ returns the sites coordinates that minimise the
time-average of $\bm{\mathcal{C}}_{ij}$, $\Omega_{ij}$, without any extra
constraints. The second site we pick does not, however, have the same degrees
of freedom, and we want to minimise both $\Omega_{ij}$ and the site-pair
correlation with $\mathcal{S}_{1}$. Hence, we define a new, quadratic
(positive, semi-definite) function, $g$, as a function of spatial coordinates
$(x_{i},y_{j})$ such that

$\displaystyle
g_{ij}=\Big{\|}\omega_{0}\Omega_{ij}+\omega_{1}r_{ij}(\bm{\mathcal{C}}_{1},\bm{\mathcal{C}}_{ij})\Big{\|}^{2},$
(12)

where $\omega_{0}$ and $\omega_{1}$ are weights that determine what we deem
more important – the cloud cover between sites, or maintaining an uncorrelated
network, and $\|a+b\|^{2}=a^{2}+b^{2}$ is the two-norm squared. Note that
$g_{ij}=0$ at its global minima, which corresponds to a site with no average
cloud cover, and no correlation with $\mathcal{S}_{1}$, and
$g_{ij}=\omega_{1}^{2}+\omega_{2}^{2}$ at its maxima, describing a site that
is on average always cloudy, and strongly correlated (either positive or
negative888Note that we show that both types of correlations reduce the
performance of the network in Appendix C) with $\mathcal{S}_{1}$. Hence,
$\mathcal{S}_{2}$ is

$\displaystyle\mathcal{S}_{2}(x_{i},y_{j})=\text{argmin}_{ij}\left\\{g_{ij}\right\\}.$
(13)

Now we must generalise Equation 12 for a general
$\mathcal{S}_{k}(x_{i},y_{j})$, noting that we ought to be always able to
minimise $g_{ij}$ to find any next site on the network. The procedure for $N$
sites is simply,

$\displaystyle g_{ij}(N)=$
$\displaystyle\Big{\|}\omega_{0}\Omega_{ij}+\frac{1}{N}\sum_{k=1}^{N}\omega_{k}r_{ij}(\bm{\mathcal{C}}_{k},\bm{\mathcal{C}}_{ij})\Big{\|}^{2},$
(14) $\displaystyle=$
$\displaystyle\Big{\|}\omega_{0}\Omega_{ij}+\frac{1}{N}\left[\omega_{1}r_{ij}(\bm{\mathcal{C}}_{1},\mathcal{C}_{ij})\right.+$
$\displaystyle\left.\omega_{2}r_{ij}(\bm{\mathcal{C}}_{2},\bm{\mathcal{C}}_{ij})+\ldots+\omega_{N}r_{ij}(\bm{\mathcal{C}}_{N},\bm{\mathcal{C}}_{ij})\right]\Big{\|}^{2},$
(15)

which describes our selection criteria, $g_{ij}(N)$ as a linear combination of
time-averaged cloud cover through the first term, and pair-wise site
correlations, each with a weighting
$\left\\{\omega_{1},\ldots,\omega_{N}\right\\}$. This maintains $g_{ij}(N)=0$
as a global minima, but now
$g_{ij}(N)=\omega_{0}^{2}+(\omega_{1}^{2}+\omega_{2}^{2}+\ldots+\omega_{N}^{2})/N^{2}$
is the global maximum. Hence, to generate the first site in $\mathcal{N}$ we
use Equation 11, and then to successively generate $n\in N$ following sites we
use Equation 14, i.e.,

$\displaystyle\mathcal{S}_{n}(x_{i},y_{j})=\text{argmin}_{i,j}\begin{cases}\left\\{\Omega_{ij}\right\\},&\text{
if }n=1,\\\ \left\\{g_{ij}(n)\right\\},&\text{ if }n>1,\end{cases}$ (16)

to populate $\mathcal{N}$. Note that this simple method can be easily
generalised for including any other site characteristics that may help
determine a robust $\mathcal{N}$. This method can also be used to optimally
add OGS to an existing network, explored in Appendix A.

We initially set all $\omega=1$, however subsection 6.3 of Section 6 explores
the use of a $\omega$ function to further optimise for a given orbital
inclination. The $\omega$ terms can also be used to mask for regions-of-
interest. The result of running this method with eight years of daily cloud
data from Himawari-8 over Australia is shown in Figure 7 for $N=8$. New
Zealand is ignored for simplicity in this analysis and to reduce computational
demand given the number of correlation calculations.

Figure 7: Distribution of $g_{ij}(N=8)$ function, Equation 12, and the minimum
for linear combinations of availability and spatial correlation from
Himawari-8 cloud retrieval data. Eight optimally selected sites,
$\mathcal{N}_{N=8}$, are highlighted with white dots and indexed from $k=1$.

We also reduce the spatial resolution of Himawari-8 data for this analysis to
25 km from 5 km to further reduce computational effort. Figure 7 is dark where
$g_{ij}$ is close to $0$, and the process of minimising $g_{ij}$ selects these
areas. Each $\mathcal{S}_{k}$ also imposes a structure to the surface of
$g_{ij}$ related to the $r_{ij}(\mathcal{C}_{k},\mathcal{C}_{ij})$ functions
in Figure 4. We see from Figure 7 that sequentially minimising $g_{ij}(N)$
picks locations for $\mathcal{N}_{N=8}$ which are well-spaced and in areas of
low cloud cover, tending to be in extremely arid parts of Australia.

The effectiveness of our optimisation technique can then be judged by
comparing with the extended AOGSN, a more realistic example network. Table 2
uses the mean of $A$ across each of the sites, $\bar{A}$, $\bar{|r|}$
(Equation 3), and $p_{\text{avail}}(M=0)$ to compare the AOGSN with
$\mathcal{N}_{N=8}$. Table 2 shows that for all metrics $\mathcal{N}_{N=8}$
outperforms the extended AOGSN. We also see in Table 2 that
$\mathcal{N}_{N=8}$ has an order of magnitude lower chance of a network-wide
outage. However, the metrics we employ in Table 2 do not properly represent
network capability upon comparison, e.g. outage probability for both
$\mathcal{N}_{N=8}$ and the extended AOGSN seem sufficient for a reliable
network. In the next sections we use orbital simulations to provide a more
detailed view of network capability.

Table 2: Comparison between base and full AOGSN with optimised network,
$\mathcal{N}$.

Network |  $\bar{A}$ |  $\bar{|r|}$ |  $p_{\text{avail}}(M=0)$
---|---|---|---
AOGSN Extended |  $0.69\pm 0.09$ |  $0.08\pm 0.11$ |  $2.4\times 10^{-4}$
AOGSN Base |  $0.67\pm 0.06$ |  $0.073\pm 0.076$ |  $0.064$
$\mathcal{N}_{N=8}$ |  $0.83\pm 0.04$ |  $0.069\pm 0.060$ |  $1.47\times 10^{-5}$

  * •

Notes. $\bar{A}$ is statistical mean of $A$ across all sites in a given
network and $\bar{|r|}$ is determined from Equation 3 shown with $\pm 1\sigma$
from both calculations. $p_{\text{avail}}(M=0)$ is from the Monte Carlo method
for sampling the availability PDF, described in Section 4 The uncertainties of
$A$, as with Fig 3, represent seasonal variation and satellite-to-satellite
variability. Uncertainties of $\bar{|r|}$ is the $\sigma$ from all values in
the matrix used to compute $\bar{|r|}$ described in Equation 3.

##  6 LEO satellite link modelling

We conduct orbital propagation simulations that can be combined with
availability measurements to estimate the data capacity for different network
configurations. We do this by using the a.i.-solutions FreeFlyer software to
simulate LEO satellite links to the AOGSN. Coverage as a function of orbital
height ($h_{\text{orbit}}$) and minimum OGS elevation angle
($\theta_{\text{elevation}}$) can be used to estimate the number and duration
of passes a satellite has with the network. While GEO or deep space
applications are likely to maintain mutual visibility across the entire
network, LEO satellites will typically only see one node at a time, so an
increase in number of passes directly corresponds to an increased number of
links and average data capacity.

Figure 8 shows coverage for ground stations in the AOGSN with
$\theta_{\text{elevation}}\geq 30^{\circ}$ and $\theta_{\text{elevation}}\geq
45^{\circ}$ to a satellite with $h_{\text{orbit}}=500$ km.
$h_{\text{orbit}}=500$ km is selected as a typical orbit for LEO satellites.
The area covered by the network in Figure 8 is large, highlighting how an
Australasian network can provide ground segment support over a significant
area, continuing a well-established role in radio ground segment operations
[37]. However, the overlapping circles in Figure 8 represent areas of mutual
coverage, implying the proposed extended AOGSN has some redundancy not present
in the optimised network, $\mathcal{N}_{N=8}$.

Figure 8: Coverage for the ground stations in the extended AOGSN with
$\theta_{\text{elevation}}>30^{\circ}$ (blue) and
$\theta_{\text{elevation}}>45^{\circ}$ (green) to satellite with
$h_{\text{orbit}}=500$ km.

###  6.1 Link duration analysis

We conduct further orbital simulations to determine average link durations
$\tau(i)$ between LEO orbits and the extended AOGSN, i.e., the sum of passes
and their duration. The parameter space of orbital inclination $i$ is modelled
by scaling a fiducial ISS orbit to $530$ km (common height for LEO
communication satellites), and then orientating it from
$i=20^{\circ}-100^{\circ}$. These orbits are propagated for one year to create
a robust statistical sample. This method may overestimate $\tau$ because it
includes short passes that glance $\theta_{elevation}=30^{\circ}$ which would
typically be ignored, and passes close to the Sun were not discarded. A single
time step of this simulation is visible in Figure 9, where sample orbits are
shown in blue, and site coverage at $i=30^{\circ}$ in green. For each site in
the AOGSN, all passes and their durations are totalled to determine $\tau(i)$
for $530$ km orbital height and a minimum $\theta_{elevation}=30^{\circ}$.

Figure 9: A single time-step from a.i.-solutions FreeFlyer orbital simulations
of AOGSN and satellites with inclination $i=10^{\circ}-100^{\circ}$ and
$h_{\text{orbit}}=530$ km. Blue traces are recent satellite orbital
trajectories. Green circles correspond to coverage for
$\theta_{\text{elevation}}\geq 30^{\circ}$. Note that the Dongara site is not
visible in this sample.

Figure 10 shows the per day average of $\tau(i)$ from the year-long orbital
simulation results. Each site has a similar curve, with a peak around a sites
respective latitude, $L$. This relationship $\max\\{\tau(i)\\}\propto L$
demonstrates that networks and OGS placement can be optimised for a particular
orbit or a given constellation, e.g. southern latitude (for the southern
hemisphere) sites are preferable to link with polar orbiting satellites.

Figure 10: The daily average link duration, $\tau(i)$ (number of seconds/day)
between each ground station (different colours) of the extended AOGSN and
satellites with orbital inclinations from $20^{\circ}-100^{\circ}$.
Availability via the cloud cover (see Section 2) is not considered for these
estimates. Orbital parameters and results are based on a fiducial ISS orbit
scaled to $530$ km and propagated for one year.

###  6.2 Latitude correction for network optimisation

Figure 10 shows the strong relation between site latitude, $L$, and $\tau(i)$.
A network optimised for availability and diversity preferences the arid
northern portion of Australia, such as $\mathcal{N}_{N=8}$ from Section 5,
which may not be ideal for maximising $\tau$ with a given satellite or
constellation. We therefore adjust the method from Section 5 to produce an
optimised network for a specific $i=90^{\circ}$ polar orbit satellite. Many
current FSOC satellites are approximately in polar orbits, and Earth
observation satellites, which may one day be supported by an FSOC ground
segment for improved downlinks, typically use a Sun-synchronous near-polar
orbit.

To optimise for an $i=90^{\circ}$ orbit we look at the $\tau(i=90^{\circ})$
values for each AOGSN site in Figure 10, and deduce a relation empirically
from our simulations. If $L$ is far from $i$, we can see that $\tau(i)$ is
asymptotic, which is satisfied for all sites in the extended AOGSN. We
therefore record $\tau(i=90^{\circ})$ as a function of $L$ and normalise the
results,

$\displaystyle\omega(L)=0.00745L+1$ (17)

to create an expression that describes $\tau(i=90^{\circ})\propto L$.

$\omega(L)$ is normalised to the equator ($L=0$) with a minimum of $0.33$ at
the south pole ($L=-90$). We can then use $\omega(L)$ as a weighting function
for all $\omega_{k}$ in Equation 12, which adds a linear North-South gradient
for preferencing southern sites to the surface, $g_{ij}(N)$, being minimised
by the algorithm in Section 5. We re-run the optimisation to generate
$\mathcal{M}_{N=8}$, the set of sites selected after applying this latitude-
correction which is shown in Figure 11. We see very little difference between
Figure 11 ($\mathcal{M}_{N=8}$) and Figure 7 ($\mathcal{N}_{N=8}$) for $S_{1}$
to $S_{6}$, but with some southern shift for $S_{7}$ and $S_{8}$ due to the
new correction term, $\omega(L)$. This similarity is likely because the first
sites have such high availability that they remain optimal despite the orbital
dynamics favouring southern locations. The method we have demonstrated could
be generalised to any orbital parameters, not just $i$, and for an entire
satellite constellation.

Figure 11: Optimised selection of eight sites, $\mathcal{M}_{N=8}$, using
method from Section 5 with linear latitude weighting function, Equation 17, to
optimise specifically for $i=90^{\circ}$ orbits.

###  6.3 Network capacity

We finally combine work on availability, $\tau(i)$ modelling, and optimisation
techniques by constructing a general network capacity for LEO satellites. The
availability-corrected, network-wide link duration, $T(i)$, for a network of
$N$ sites is,

$T(i)=\sum_{k=1}^{N}A_{k}\tau_{k}(i).$ (18)

We calculate $T(i)$ for four networks: (1) the base AOGSN, (2) extended AOGSN,
(3) $\mathcal{N}_{N=8}$, and (4) $\mathcal{M}_{N=8}$. Figure 12 shows $T(i)$
for these four networks, and an estimate of the daily data capacity with a $5$
Gbps bit rate, possible with future on-off keying satellite terminals [38,
39]. Note that this expression for $T(i)$ overestimates total link time by
including scenarios where a LEO can maintain mutual visibility with $>1$
available sites, by including both links. In practice, the satellite must
choose which OGS to link with, and the chance of both being available
simultaneously is related to $r$ for that pair, e.g. from Figure 5. Estimating
network capacity in this manner does not consider potential outages or
degradation of bit-rate due to turbulence, so it is only a guide of system
performance, for a given bit-rate.

Figure 12 shows that both the AOGSN and algorithm-optimised networks can
feasibly support many terabits of data transfer per day with a LEO satellite
and current capability [39, 38]. An general improvement can be seen in Figure
12 between the extended AOGSN and optimised networks, $\mathcal{N}_{N-8}$ and
$\mathcal{M}_{N=8}$. $\mathcal{N}_{N=8}$ and $\mathcal{M}_{N=8}$ have a $13\%$
and $16\%$ increased data capacity compared with the extended AOGSN
respectively for $i=100^{\circ}$, approximately a Sun-synchronous polar orbit.
Therefore, the minor difference of $S_{7},S_{8}$ $\in\mathcal{M}_{N=8}$,
optimised with the $\omega(L)$ weighting function, shows only a marginal
improvement of $3\%$ over $\mathcal{N}_{N=8}$ for highly inclined orbits.

We may be interested in the case of network capacity for a LEO constellation
uniformly distributed in $i$. Integrating $T(i)$ over $i$ for all networks in
Figure 12 allows to compare this general case. Upgrading from the $N=3$ base
AOGSN to $N=8$ extended AOGSN corresponds to a $293\%$ increase of data
capacity when integrating $T(i)$ for the uniform case. This also shows a
$27\%$ and $20\%$ efficiency gain for $\mathcal{N}_{N=8}$ and
$\mathcal{M}_{N=8}$ respectively, compared to the extended AOGSN. The
optimisation from Section 5 therefore yields a network that is superior in
data capacity to the more realistic AOGSN. We could also use this method to
optimally select additional sites for a network and estimate the increased
capacity from doing so.

Figure 12: Network link duration, $T(i)$ (Equation 18), and data capacity
estimates for the AOGSN (extended, solid; base, dashed) and two networks
optimised for maximising availability and diversity. Optimised networks
$\mathcal{M}_{N=8}$ (dashed-dotted, subsection 6.2) and $\mathcal{N}_{N=8}$
(dotted with triangles, Section 5) are with and without latitude-correction
for polar orbits respectively. Data estimation uses $5$ Gbps as the fiducial
bitrate, multiplied by $T(i)$ for current on-off keying capability, assuming a
uniform atmospheric channel for every site. Figure 13: Number of visible AOGSN
(extended, solid line; base, dashed line) nodes (blue) and corresponding
outage probability (green), $p_{\text{avail}}(M=0)$, for GEO satellites as a
function of longitude. GEO visibility to AOGSN ground stations (orange) is
with $\theta_{\text{elevation}}>30^{\circ}$. Outage probabilities are
estimated by sampling the joint probability distribution of site availability,
as in Section 4.

##  7 GEO visibility analysis

While we have extensively studied the case of how the AOGSN can support LEO
satellites, a large and diverse optical network can provide very high
reliability to GEO or deep space if there is significant mutual visibility to
numerous OGS. GEO optical feeder links to ground may become an important
future application of FSOC networks, where ground segment reliability is
critical [40, 41]. The Asia-Pacific region may soon host a GEO satellite with
an FSOC terminal to demonstrate this capability, which an Australasian network
could reliably support [42, 39].

We therefore use orbital simulations with a.i.-solutions’ FreeFlyer to
determine GEO visibility to each ground station in the AOGSN with
$\theta_{elevation}>30^{\circ}$. For simplicity, GEO orbital placement is
considered only a function of longitude with $h_{\text{orbit}}=35,786$ km and
$i\lesssim 1^{\circ}$. Hence, GEO longitude corresponds to visibility with
$N\leq 8$ nodes, allowing estimates of $p_{\text{avail}}(M=0)$, described in
Section 4, as a function of longitude. Figure 13 shows the number of nodes in
the extended AOGSN visible to a GEO as a function of longitude. A portion of
the AOGSN is therefore visible to a GEO between different longitude bounds.
For each segment, e.g. $N=8$ from $140^{\circ}$ to $160^{\circ}$, Figure 13
shows $p_{\text{avail}}(M=0)$ for that set of sites. The extended AOGSN is
three times longer in longitude than latitude, which helps maintain
$p_{\text{avail}}(M=0)<0.01$ over $\approx 80^{\circ}$ of longitude as shown
in Figure 13. This example network could then feasibly provide reliable
support to GEO satellites over nearly $1/4$ of the planet. The low number of
visible nodes to the East and West of Australia in Figure 13 indicates
extensions in these directions are ideal if minimising GEO communication
outages are a high priority. Placing an OGS in the Polynesian region, e.g.
Tahiti, could be an ideal extension by allowing visibility for North American
GEO satellites such as the Laser Communications Relay Demonstration [43].
Figure 13 further highlights the improvements between the base and extended
AOGSN.

##  8 Conclusion

###  8.1 Summary and key results

We have developed and shown use cases for a number of novel methods for
assessing and optimising FSOC ground segment networks, focusing, but not
limited to, the availability and site-wise correlations within the
Australasian region. An example network called the AOGSN is proposed as a
realistic representation of an Australasian network, with base ($N=3$) and
extended ($N=8$) variants. We furthermore perform detailed case studies on the
base and extended AOGSN, highly-relevant to the future of optical
communications within the Australasian region, comparing it to theoretically
optimised networks throughout the study.

For determining availability statistics of the networks, we compare between a
range of remote sensing data from five Earth observation satellite sources,
but particularly focus on the highly spatially and temporally resolved cloud
cover data from the Himawari-8 GEO satellite. We model the cloud cover as a
Bernoulli process, which provides us a mathematical formalism for the site
availability, and joint availability probability distribution function, which
we construct and sample to access outage probabilities. We provide a detailed
analysis of the spatially resolved correlations between pair-wise ground
stations, and between ground stations and arbitrary sites. These statistics
are combined with orbital simulations to estimate link duration and network
capacity statistics. We list the key results below:

  * •

We estimate a $69\%$ average site availability, and an average correlation of
$\bar{|r|}=0.088$, for the extended AOGSN (comparable correlations to an
ideal, optimised network; see Table 2), demonstrating the large (spatially
uncorrelated) and arid nature of the Australasian region.

  * •

Using site availability and correlation measurements and the method outlined
in [17], in Section 4 we construct the network availability probability
distribution function for an arbitrary sized Australasian network, with
arbitrary correlations between sites. The analysis of this correlated network
shows the extended AOGSN has $p_{\text{avail}}(M=0)=2.4\times 10^{-4}$, i.e.
at any given time there is close to a one-in-five thousand chance of suffering
a network-wide outage (shown in Figure 6, alongside statistics for the AOGSN
base and ideal, optimised networks), hence a large and diversified network
such as the AOGSN is reliable enough to overcome the reliability challenge for
FSOC space-to-ground communication.

  * •

We develop a spatially-resolved method for optimising the positions of ground
stations based on the network availability and diversity, without site
preselection. This method can be used to develop new or extend existing
networks, which we show in Section 5 and Appendix A, respectively, and to
generate an optimal network to provide an ideal reference for commissioned
networks (e.g., for comparing with the preselected AOGSN).

  * •

In Section 6 we combine orbital LEO simulations with previous our availability
models to estimate network capacity as a function of orbital inclination,
showing that the extended AOGSN network is close to optimal (compared to our
ideal, optimised network) and both the AOGSN and optimal networks are capable
of supporting bitrates of $\gtrsim 10$ Tb/day to a LEO.

###  8.2 Limitations of our study and future work

In this study we consider only a simple availability model based on time-
averaged cloud cover. Clearly, as shown in Figure 3, there is seasonal
variation (which can lead up to a factor of $\approx 2$ variation in the
availability) and temporal correlation that we do not account for in our
study. We will leave a detailed analysis of how the performance of FSOC
networks in the Australasian network will change with seasonal, climatic and
temporal correlations for future work. Cirrus clouds, and other thin icy
clouds, could be more rigorously discarded from the total cloud fraction, as
they usually do not create an outage. A method for differentiating the liquid
cloud fraction, i.e. approximately thick clouds, from the total cloud fraction
is outlined in [44] and could be adapted in future studies that want to add
additionally complexity to the models we have presented.

Furthermore, the atmospheric channel will also constrain availability, which
has been explored theoretically in other network diversity analyses [11, 3].
Very arid and flat sites with low cloud cover likely have strong and
detrimental turbulence conditions, hence, the results presented in our study
represent an upper bound of the site availabilities. This consideration may be
especially important in the context of sensitive coherent protocols, which may
require adaptive optics [45]. Network degradation due to turbulence will be
highly system specific, so estimating system performance should be done on a
case-by-case basis. We plan to deploy atmospheric turbulence instruments in
prospective sites as part of future investigations to explore OGS placement
for the Australasian region, with a focus on the suitability of low cloud
cover desert and arid sites. However, for a network model that retains spatial
freedom and without an experimental campaign at selected sites, there is no
feasible or reasonable way to include turbulence in the present, spatially
resolved outage model. Well-validated spatial estimates of daytime turbulence
do not exist and testing emerging models (e.g. [46, 47]), will require
extensive experimental data for future models that include these effects.

Our optimisation method described in Section 5 relies upon only two site
characteristics: (1) the availability and (2) the spatially-resolved cloud
cover correlation, both of which are equally weighted and without any spatial
(latitude, longitude or altitude) constraints. The latitude-correction case
studied in subsection 6.2 of Section 6 demonstrates the flexibility of this
method by optimising for a specific orbit. The theoretical outage model
involving the atmospheric channel developed in [3] highlights that our cloud-
only model, and proceeding literature, can only provide an upper bound of
availability by not including turbulence-induced outages [10].

## Funding

M. B. acknowledges funding from the Australian Government Research Training
Program and the CSIRO postgraduate support scholarship. J. R. B. acknowledges
financial support from the Australian National University, via the Deakin Ph.D
and Dean’s Higher Degree Research (theoretical physics) Scholarships, the
Australian Government via the Australian Government Research Training Program
and the Australian Capital Territory Government funded Fulbright scholarship.
N. J. R.and J. E. C are supported by a New Zealand-DLR Joint Research
Programme funded by New Zealand’s Ministry of Business, Innovation and
Employment.

## Acknowledgements

The Terra/MODIS, Aqua/MODIS, Suomi NPP/VIIRS, and Aqua/AIRS datasets were
acquired from the Level-1 and Atmosphere Archive & Distribution System (LAADS)
Distributed Active Archive Center (DAAC), located in the Goddard Space Flight
Center in Greenbelt, Maryland. Himawari-8 data used for this paper has been
provided by the Earth Observation Research Centre of the Japan Aerospace
Exploration Agency. Data from the Modern-Era Retrospective analysis for
Research and Applications, Version 2 (MERRA-2) was provided by the National
Aeronautics and Space Administration Global Modelling and Assimilation Office,
located in the Goddard Space Flight Center in Greenbelt, Maryland. We also
thank the anonymous reviewers who helped enhance the clarity and presentation
of our study.

We acknowledge the open-source data analysis and visualisation software used
in this study: numpy [48, 49], matplotlib [50], scipy [35], cartopy [51].

## References

  * [1] F. Bennet, K. Ferguson, K. Grant, E. Kruzins, N. Rattenbury, and S. Schediwy, “An Australia/New Zealand optical communications ground station network for next generation satellite communications,” in _Free-Space Laser Communications XXXII,_ vol. 11272 (SPIE, 2020), p. 1127202.
  * [2] S. Piazzolla and S. Slobin, “Statistics of link blockage due to cloud cover for free-space optical communications using NCDC surface weather observation data,” in _Free-Space Laser Communication Technologies XIV,_ vol. 4635 (SPIE, 2002), pp. 138–149.
  * [3] E. Erdogan, I. Altunbas, G. K. Kurt, M. Bellemare, G. Lamontagne, and H. Yanikomeroglu, “Site diversity in downlink optical satellite networks through ground station selection,” IEEE Access 9, 31179–31190 (2021).
  * [4] I. del Portillo, M. Sanchez, B. Cameron, and E. Crawley, “Optimal location of optical ground stations to serve LEO spacecraft,” in _Proc. IEEE Aerosp. Conf._ , (2017), pp. 1–16.
  * [5] M. S. Net, I. Del Portillo, E. Crawley, and B. Cameron, “Approximation methods for estimating the availability of optical ground networks,” Journal of Optical Communications and Networking 8, 800–812 (2016).
  * [6] N. K. Lyras, C. N. Efrem, C. I. Kourogiorgas, and A. D. Panagopoulos, “Medium earth orbit optical satellite communication networks: Ground terminals selection optimization based on the cloud-free line-of-sight statistics,” International Journal of Satellite Communications and Networking 37, 370–384 (2019).
  * [7] S. Walsh, A. Frost, W. Anderson, T. Digney, B. Dix-Matthews, D. Gozzard, C. Gravestock, L. Howard, S. Karpathakis, A. McCann _et al._ , “The Western Australian optical ground station,” 38th International Communications Satellite Systems Conference (ICSSC 2021) (2021).
  * [8] M. Birch, N. Martínez, F. Bennet, M. Copeland, and D. Grosse, “The Mount Stromlo Optical Communication Ground Station,” in _2022 IEEE International Conference on Space Optical Systems and Applications (ICSOS),_ (IEEE, 2022), pp. 134–137.
  * [9] K. Mudge, B. Clare, E. Jager, V. Devrelis, F. Bennet, M. Copeland, N. Herrald, I. Price, G. Lechner, J. Kodithuwakkuge _et al._ , “DSTG Laser Satellite Communications-Current Activities and Future Outlook,” in _2022 IEEE International Conference on Space Optical Systems and Applications (ICSOS),_ (IEEE, 2022), pp. 17–21.
  * [10] C. Fuchs and F. Moll, “Ground station network optimization for space-to-ground optical communication links,” Journal of Optical Communications and Networking 7, 1148–1159 (2015).
  * [11] N. K. Lyras, C. N. Efrem, C. I. Kourogiorgas, and A. D. Panagopoulos, “Optimum monthly based selection of ground stations for optical satellite networks,” IEEE Communications Letters 22, 1192–1195 (2018).
  * [12] N. Martínez, L. F. Rodríguez-Ramos, and Z. Sodnik, “Toward the uplink correction: application of adaptive optics techniques on free-space optical communications through the atmosphere,” Optical Engineering 57, 076106–076106 (2018).
  * [13] J. Osborn, M. J. Townson, O. J. Farley, A. Reeves, and R. M. Calvo, “Adaptive optics pre-compensated laser uplink to leo and geo,” Optics Express 29, 6113–6132 (2021).
  * [14] C. M. Schieler, K. M. Riesing, B. C. Bilyeu, J. S. Chang, A. S. Garg, N. J. Gilbert, A. J. Horvath, R. S. Reeve, B. S. Robinson, J. P. Wang _et al._ , “On-orbit demonstration of 200-gbps laser communication downlink from the tbird cubesat,” in _Free-Space Laser Communications XXXV,_ vol. 12413 (SPIE, 2023), p. 1241302.
  * [15] L. C. Roberts, S. R. Meeker, J. Tesch, J. C. Shelton, J. E. Roberts, S. F. Fregoso, T. Troung, M. Peng, K. Matthews, H. Herzog _et al._ , “Performance of the adaptive optics system for laser communications relay demonstration’s ground station 1,” Applied Optics 62, G26–G36 (2023).
  * [16] K. Bessho, K. Date, M. Hayashi, A. Ikeda, T. Imai, H. Inoue, Y. Kumagai, T. Miyakawa, H. Murata, T. Ohno _et al._ , “An introduction to Himawari-8/9—Japan’s new-generation geostationary meteorological satellites,” Journal of the Meteorological Society of Japan. Ser. II 94, 151–183 (2016).
  * [17] J. H. Macke, P. Berens, A. S. Ecker, A. S. Tolias, and M. Bethge, “Generating Spike Trains with Specified Correlation Coefficients,” Neural Comput. 21, 397–423 (2009).
  * [18] S. Poulenard, M. Crosnier, and A. Rissons, “Ground segment design for broadband geostationary satellite with optical feeder link,” Journal of Optical Communications and Networking 7, 325–336 (2015).
  * [19] R. A. Venkat and D. W. Young, “Cloud-free line-of-sight estimation for free space optical communications,” in _Atmospheric Propagation X,_ vol. 8732 (SPIE, 2013), pp. 17–26.
  * [20] S. Platnick, K. G. Meyer, M. D. King, G. Wind, N. Amarasinghe, B. Marchant, G. T. Arnold, Z. Zhang, P. A. Hubanks, R. E. Holz _et al._ , “The MODIS cloud optical and microphysical products: Collection 6 updates and examples from Terra and Aqua,” IEEE Transactions on Geoscience and Remote Sensing 55, 502–525 (2016).
  * [21] H. H. Aumann, M. T. Chahine, C. Gautier, M. D. Goldberg, E. Kalnay, L. M. McMillin, H. Revercomb, P. W. Rosenkranz, W. L. Smith, D. H. Staelin _et al._ , “AIRS/AMSU/HSB on the Aqua mission: Design, science objectives, data products, and processing systems,” IEEE Transactions on Geoscience and Remote Sensing 41, 253–264 (2003).
  * [22] S. Platnick, K. Meyer, G. Wind, R. E. Holz, N. Amarasinghe, P. A. Hubanks, B. Marchant, S. Dutcher, and P. Veglio, “The NASA MODIS-VIIRS continuity cloud optical properties products,” Remote sensing 13, 2 (2020).
  * [23] R. A. Frey, S. A. Ackerman, R. E. Holz, S. Dutcher, and Z. Griffith, “The continuity MODIS-VIIRS cloud mask,” Remote Sensing 12, 3334 (2020).
  * [24] R. Gelaro, W. McCarty, M. J. Suárez, R. Todling, A. Molod, L. Takacs, C. A. Randles, A. Darmenov, M. G. Bosilovich, R. Reichle _et al._ , “The modern-era retrospective analysis for research and applications, version 2 (MERRA-2),” Journal of climate 30, 5419–5454 (2017).
  * [25] M. H. England, C. C. Ummenhofer, and A. Santoso, “Interannual rainfall extremes over southwest Western Australia linked to Indian Ocean climate variability,” Journal of Climate 19, 1948–1969 (2006).
  * [26] J. R. Norris, R. J. Allen, A. T. Evan, M. D. Zelinka, C. W. O’Dell, and S. A. Klein, “Evidence for climate change in the satellite cloud record,” Nature 536, 72–75 (2016).
  * [27] G. Tselioudis, B. R. Lipat, D. Konsta, K. M. Grise, and L. M. Polvani, “Midlatitude cloud shifts, their primary link to the Hadley cell, and their diverse radiative effects,” Geophysical Research Letters 43, 4594–4601 (2016).
  * [28] J. Choi, S.-W. Son, J. Lu, and S.-K. Min, “Further observational evidence of Hadley cell widening in the Southern Hemisphere,” Geophysical Research Letters 41, 2590–2597 (2014).
  * [29] H. Nguyen, C. Lucas, A. Evans, B. Timbal, and L. Hanson, “Expansion of the Southern Hemisphere Hadley cell in response to greenhouse gas forcing,” Journal of Climate 28, 8067–8077 (2015).
  * [30] H. Ishida and T. Y. Nakajima, “Development of an unbiased cloud detection algorithm for a spaceborne multispectral imager,” Journal of Geophysical Research: Atmospheres 114 (2009).
  * [31] H. Iwabuchi, N. S. Putri, M. Saito, Y. Tokoro, M. Sekiguchi, P. Yang, and B. A. Baum, “Cloud property retrieval from multiband infrared measurements by Himawari-8,” Journal of the Meteorological Society of Japan. Ser. II (2018).
  * [32] M. S. Awan, E. Leitgeb, B. Hillbrand, F. Nadeem, M. Khan _et al._ , “Cloud attenuations for free-space optical links,” in _2009 International Workshop on Satellite and Space Communications,_ (IEEE, 2009), pp. 274–278.
  * [33] D. Tjøstheim, H. Otneim, and B. Støve, “Statistical dependence: Beyond Pearson’s $\rho$,” arXiv e-prints arXiv:1809.10455 (2018).
  * [34] B. Dai, S. Ding, and G. Wahba, “Multivariate bernoulli distribution,” Bernoulli 19, 1465–1483 (2013).
  * [35] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright _et al._ , “SciPy 1.0: fundamental algorithms for scientific computing in Python,” Nature methods 17, 261–272 (2020).
  * [36] M. Toyoshima, Y. Munemasa, H. Takenaka, Y. Takayama, Y. Koyama, H. Kunimori, T. Kubooka, K. Suzuki, S. Yamamoto, S. Taira _et al._ , “Introduction of a terrestrial free-space optical communications network facility: IN-orbit and Networked Optical ground stations experimental Verification Advanced testbed (INNOVA),” in _Free-Space Laser Communication and Atmospheric Propagation XXVI,_ vol. 8971 (SPIE, 2014), pp. 207–212.
  * [37] I. del Portillo, M. Sanchez, B. Cameron, and E. Crawley, “Architecting the ground segment of an optical space communication network,” in _2016 IEEE Aerospace Conference,_ (IEEE, 2016), pp. 1–13.
  * [38] C. Schmidt, B. Rödiger, J. Rosano, C. Papadopoulos, M.-T. Hahn, F. Moll, and C. Fuchs, “DLR’s Optical Communication Terminals for CubeSats,” in _2022 IEEE International Conference on Space Optical Systems and Applications (ICSOS),_ (IEEE, 2022), pp. 175–180.
  * [39] A. Carrasco-Casado, K. Shiratama, P. V. Trinh, D. Kolev, Y. Munemasa, T. Fuse, H. Tsuji, and M. Toyoshima, “Development of a miniaturized laser-communication terminal for small satellites,” Acta Astronautica (2022).
  * [40] A. M. Bonnefois, J.-M. Conan, C. Petit, C. Lim, V. Michau, S. Meimon, P. Perrault, F. Mendez, B. Fleury, J. Montri _et al._ , “Adaptive optics pre-compensation for GEO feeder links: the FEEDELIO experiment,” in _International Conference on Space Optics—ICSO 2018,_ vol. 11180 (SPIE, 2019), pp. 889–896.
  * [41] N. Védrenne, A. Montmerle-Bonnefois, C. Lim, C. Petit, J.-F. Sauvage, S. Meimon, P. Perrault, F. Mendez, B. Fleury, J. Montri _et al._ , “First experimental demonstration of adaptive optics pre-compensation for GEO feeder links in a relevant environment,” in _2019 IEEE International Conference on Space Optical Systems and Applications (ICSOS),_ (IEEE, 2019), pp. 1–5.
  * [42] D. Kolev, K. Shiratama, A. Carrasco-Casado, Y. Saito, Y. Munemasa, J. Nakazono, P. V. Trinh, H. Kotake, H. Kunimori, T. Kubooka _et al._ , “Status Update on Laser Communication Activities in NICT,” in _2022 IEEE International Conference on Space Optical Systems and Applications (ICSOS),_ (IEEE, 2022), pp. 36–39.
  * [43] D. J. Israel, B. L. Edwards, and J. W. Staren, “Laser Communications Relay Demonstration (LCRD) update and the path towards optical relay operations,” in _2017 IEEE Aerospace Conference,_ (IEEE, 2017), pp. 1–6.
  * [44] C. Listowski, J. Delanoë, A. Kirchgaessner, T. Lachlan-Cope, and J. King, “<? xmltex 3mm?> antarctic clouds, supercooled liquid water and mixed phase, investigated with dardar: geographical and seasonal variations,” Atmospheric Chemistry and Physics 19, 6771–6808 (2019).
  * [45] M. Gregory, D. Troendle, G. Muehlnikel, F. Heine, R. Meyer, M. Lutzer, and R. Czichy, “Three years coherent space to ground links: performance results and outlook for the optical ground station equipped with adaptive optics,” in _Free-Space Laser Communication and Atmospheric Propagation XXV,_ vol. 8610 (SPIE, 2013), pp. 17–29.
  * [46] S. Basu, J. Osborn, P. He, and A. DeMarco, “Mesoscale modelling of optical turbulence in the atmosphere: the need for ultrahigh vertical grid resolution,” Monthly Notices of the Royal Astronomical Society 497, 2302–2308 (2020).
  * [47] J. Osborn and M. Sarazin, “Atmospheric turbulence forecasting with a general circulation model for cerro paranal,” Monthly Notices of the Royal Astronomical Society 480, 1278–1299 (2018).
  * [48] T. Oliphant, “NumPy: A guide to NumPy,” USA: Trelgol Publishing (2006). [Online; accessed <today>].
  * [49] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, “Array programming with numpy,” Nature 585, 357–362 (2020).
  * [50] J. D. Hunter, “Matplotlib: A 2d graphics environment,” Computing in Science & Engineering 9, 90–95 (2007).
  * [51] Met Office, _Cartopy: a cartographic python library with a Matplotlib interface_ , Exeter, Devon (2010 - 2015).

## Appendix

##  Appendix A Optimisation for additional sites and regions

Our optimisation method described in Section 5 can be used to select
additional sites for an existing network. This is done by first constructing
the $g_{ij}(N)$ function for the initial ground stations and then adding sites
following the process used to construct $\mathcal{N}_{N=8}$ and
$\mathcal{M}_{N=8}$. We test this using the AOGSN base network as the initial
condition, and run the algorithm until it reaches $N=8$, as shown in Figure
14. This method could be particularly useful when combined with the capacity
estimation from Section 6, as the improvement from optimally adding a new OGS
to an existing network can be quantified. This application could also be
useful with additional ROI masking functions to guide site selection for
extending a given network, such as the AOGSN. Figure 14 indicates that for
optimisation, the best site to extend the base AOGSN with would be located in
northern Western Australia.

Figure 14: The same as Figure 7 but for $N=5$ with the base AOGSN ground
stations as initial conditions. Daily AHI/Himawari-8 images from 2015-2022 are
used for determining the combined $g_{ij}(N=5)+g_{\text{AOGSN}}(N=3)$
optimisation function. Figure 15: Availability as a function of month for
Mildura AOGSN site, showing the seasonal variability and each separate cloud
retrieval product (different colours) used in this paper, which are averaged
over in Figure 3.

##  Appendix B Remote sensing satellite variability

In Section 2 we outlined the method of using all five satellites from Table 1
to estimate $A$ (or $\Omega$) in a robust fashion. Figure 15 shows the
different seasonal $A$ curves from all five data products. We can see clearly
from Figure 15 that there is reasonable agreement between all satellites
despite different methods and resolution. However, Alice Springs suffered
issues with Himawari and MERRA-2, while Mt John had issues with Aqua/AIRS and
MERRA-2 due to snow and ice reflectance. In these cases, those data streams
were discarded as outliers. From Figure 15, we can see that higher resolution
data, Himawari-8, finds agreement with low resolution data, indicating that
spatial averaging is often appropriate. The reliability of spatial averaging
owes to Australia’s relatively homogeneous terrain, and could cause errors in
other environments. High resolution swathes from the Aqua, Terra, and Suomi
NPP LEO satellites can be used if the extra resolution is required.

Figure 16: The same as Figure 6 but for three networks which admit to the
AOGSN extended configuration with negative correlation (solid; $r=-0.2$), no
correlation (dashed; $r=0.0$) and positive correlation (dotted; $r=0.2$)
between ground stations.

##  Appendix C A network with negative correlations between ground stations

The most prominent negative correlations we observe in Figure 5 are due to
seasonal inversions. Negative correlations between network nodes intuitively
boosts diversity, i.e. in the $N=2$ case the minima of their availability
curves are out of phase. We test how negative correlations influence the
network statistics by constructing cumulative distributions functions, as in
Section 4 for a negatively correlated network $(r=-0.2)$, an uncorrelated
network $(r=0.0)$, and positively correlated network $(r=0.2)$, using the
AOGSN extended ground station configuration. We show the results in Figure 16,
where the $r=-0.2$ is shown with the solid line, $r=0.0$ with the dashed, and
$r=0.2$ with the dotted. We find that all networks have the same behaviour for
large $M$, as we showed previously in Figure 6, and discussed in Section 4.
However, for $M=0$ to $M\leq 4$, the networks show different availability
statistics. For this range of $M$, the network without correlations performs
the best, followed by the negatively correlated network and then the
positively correlated network. This shows that negatively correlated networks
admit to better availability statistics than positively correlated networks,
but networks without correlations between ground stations perform the best.