Further author information: (Send correspondence to M.S. and C.A.)
M.S.: E-mail<EMAIL_ADDRESS>Telephone: +39 049 829 3547
C.A.: E-mail<EMAIL_ADDRESS>Telephone: +33 049 829 3414

# MICADO PSF-Reconstruction work package description

Matteo Simioni  INAF Osservatorio Astronomico di Padova, Vicolo
dell’Osservatorio 5, I-35122, Padova, Italy  Carmelo Arcidiacono  INAF
Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Andrea Grazian  INAF Osservatorio Astronomico di Padova, Vicolo
dell’Osservatorio 5, I-35122, Padova, Italy  Yann Clenet  LESIA, Observatoire
de Paris, Section de Meudon 5, place Jules Janssen, F-92195 Meudon Cedex,
France  Richard Davies  MPE - Max-Planck-Institut für extraterrestrische
Physik Giessenbachstrasse 1, D-85748 Garching, Germany  Marco Gullieuszik
INAF Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Gijs Verdoes Kleijn  University of Groningen, PO Box 72, 9700
AB Groningen, the Netherlands  Fernando Pedichini  INAF - Osservatorio
Astronomico di Roma, Via Frascati 33, 00040, Monte Porzio Catone, Italy
Roland Wagner  Industrial Mathematics Institute, Johannes Kepler University
Linz, Altenberger Strasse 69, 4040 Linz, Austria  Ronny Ramlau  Industrial
Mathematics Institute, Johannes Kepler University Linz, Altenberger Strasse
69, 4040 Linz, Austria  RICAM - Johann Radon Institute for Computational and
Applied Mathematics, Altenberger Strasse 69, 4040 Linz, Austria  Werner W.
Zeilinger  Department of Astrophysics, University of Vienna,
Tuerkenschanzstrasse 17, A-1180, Wien, Austria  Fabrice Vidal  LESIA,
Observatoire de Paris, Section de Meudon 5, place Jules Janssen, F-92195
Meudon Cedex, France  Benedetta Vulcani  INAF Osservatorio Astronomico di
Padova, Vicolo dell’Osservatorio 5, I-35122, Padova, Italy  Roberto Ragazzoni
INAF Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Dipartimento di Fisica e Astronomia, Università degli studi di
Padova, via Marzolo 8, I-35131, Padova, Italy  Arnaud Sevin  LESIA,
Observatoire de Paris, Section de Meudon 5, place Jules Janssen, F-92195
Meudon Cedex, France  Bernardo Salasnich  INAF Osservatorio Astronomico di
Padova, Vicolo dell’Osservatorio 5, I-35122, Padova, Italy  Andrea Baruffolo
INAF Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Lorenzo Busoni  INAF - Osservatorio Astrofisico di Arcetri, Via
E. Fermi 5, 50125, Firenze, Italy  Simone Esposito  INAF - Osservatorio
Astrofisico di Arcetri, Via E. Fermi 5, 50125, Firenze, Italy  Éric Gendron
LESIA, Observatoire de Paris, Section de Meudon 5, place Jules Janssen,
F-92195 Meudon Cedex, France  Roberto Piazzesi  INAF - Osservatorio
Astronomico di Roma, Via Frascati 33, 00040, Monte Porzio Catone, Italy  Elisa
Portaluri  INAF - Osservatorio Astronomico d’Abruzzo, Via Mentore Maggini,
I-64100 Teramo, Italy  Anita Zanella  INAF Osservatorio Astronomico di Padova,
Vicolo dell’Osservatorio 5, I-35122, Padova, Italy  Tapio Helin  LUT
University, P.O.Box 20, FI-53851, Lappeenranta Finland  Hanindyo Kuncarayakti
Tuorla Observatory, Department of Physics and Astronomy, FI-20014 University
of Turku, Finland  Seppo Mattila  Tuorla Observatory, Department of Physics
and Astronomy, FI-20014 University of Turku, Finland  Renato Falomo  INAF
Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, I-35122,
Padova, Italy  Enrico Pinna  INAF - Osservatorio Astrofisico di Arcetri, Via
E. Fermi 5, 50125, Firenze, Italy

###### Abstract

The point spread function reconstruction (PSF-R) capability is a deliverable
of the MICADO@ESO-ELT project. The PSF-R team works on the implementation of
the instrument software devoted to reconstruct the point spread function
(PSF), independently of the science data, using adaptive optics (AO) telemetry
data, both for Single Conjugate (SCAO) and Multi-Conjugate Adaptive Optics
(MCAO) mode of the MICADO camera and spectrograph. The PSF-R application will
provide reconstructed PSFs through an archive querying system to restore the
telemetry data synchronous to each science frame that MICADO will generate.
Eventually, the PSF-R software will produce the output according to user
specifications. The PSF-R service will support the state-of-the-art scientific
analysis of the MICADO imaging and spectroscopic data.

###### keywords:

PSF Reconstruction, MICADO, MAORY, ELT, Telemetry data, WFS

##  1 INTRODUCTION

MICADO[1] is the Multi-AO Imaging Camera for Deep Observations. It is one of
the first light instruments of the ESO ELT[2, 3, 4], providing capability both
for imaging and long-slit spectroscopy at near-infrared wavelengths while
exploiting single conjugate adaptive optics[5] (SCAO) and multiconjugate
adaptive optics[6, 7, 8] (MCAO). This provides a high level of flexibility for
modes that are not offered by other ELT instruments, giving the opportunity to
reach unique and major scientific results.

The laser guide star[9, 10] MCAO system is currently being developed by the
MAORY[11] consortium, while the natural guide star (NGS) SCAO system is
developed jointly by the MICADO and MAORY consortia[12]. MICADO will interface
to the MAORY warm optical relay that re-images the telescope focus. In this
configuration, both MCAO and SCAO are available. If required for an initial
phase, MICADO will also be able to operate with just the SCAO system in a
stand-alone mode, using a simpler optical relay that interfaces directly to
the telescope. MICADO has four observing modes producing data: standard
imaging, astrometric imaging, high-contrast imaging, and spectroscopy. Each
observing mode can optionally be operated in a configuration that gathers the
data (wavefront sensor, WFS, images and telemetry) required for point spread
function reconstruction (PSF-R).

PSF-R is the process to determine the point spread function (PSF) of a science
exposure, without using the focal plane science images. The reconstructed PSF
is instead estimated using (i) images and telemetry generated by the AO
system, (ii) information by other systems monitoring the atmosphere and (iii)
telemetry from the ELT and MICADO instrument. The MICADO PSF-R is a service
that complements, but is not part of, the MICADO data processing pipeline. It
is a tool that assists data analysis rather than data reduction. But to date,
no ground-based observatory offers PSF-R as a facility tool for AO
observations. The development of the PSF-R is a collaborative effort involving
involving the ELT and MAORY as well as MICADO, since all these systems
contribute to the final measured PSF, and they all have a role in saving and
archiving the necessary data.

The MICADO PSF-R service will support standard imaging, astrometric imaging,
and spectroscopic observations. On the other hand, the PSF-R service is not
planned for high-contrast imaging observations, as these focus on subtracting
the PSF using ad-hoc techniques. The PSF-R will not support seeing enhanced
and seeing limited (without AO) observations either. Finally, PSF-R is not
suitable for non-sidereal imaging data in case a spatially-resolved moving
target is used as AO guide star for the WFS.

The accuracy of the reconstructed PSF is driven by the scientific requirements
of MICADO. The goals of the PSF-R service is to provide reconstructed PSFs
with defined metrics (Strehl ratio, encircled energy, full width at half
maximum, FWHM, ellipticity) within the reference values summarized in Table 1.
The quoted accuracies have been derived by analytic calculations and they will
be refined with simulations being carried out by the MICADO science team and
by the scientific evaluation work package (WP) of the MICADO PSF-R team.

Table 1: Error budget for the reconstructed PSF parameters. Error budget for PSF-R is indicated in the third and fourth column respectively for the nominal case and the operational case in extreme conditions (i.e. strong turbulence, faint NGS, science object far away from the science target, incomplete AO data). parameter | notes | nominal accuracy | accuracy - extreme conditions
---|---|---|---
Strehl Ratio | SCAO on-axis | $\pm 2\%$ | $\pm 10\%$
| SCAO off-axis | $\pm 10\%$ |  $\pm 20\%$ without $C_{n}^{2}$-profile
| MCAO full field | $\pm 2\%$ | $\pm 10\%$
Encircled Energy |  | $\pm 5\%$ | $\pm 10\%$
FWHM |  | $\pm 5\%$ | $\pm 10\%$
Ellipticity | PSF axis ratio at FWHM | $\pm 5\%$ | $\pm 10\%$

##  2 WORK PACKAGES AND TEAM MEMBERS

The MICADO PSF-R team works on the implementation of the instrument software
devoted to reconstruct the PSF, independently of the science data, using AO
telemetry data, both for SCAO and MCAO mode of the MICADO camera and
spectrograph. The work is distributed over $9$ WPs, listed in Table 2,
together with the institutions responsible for their development. In the
following, the goals of each WP are summarized.

Table 2: MICADO PSF-R WPs. For each WP, the responsible institution is also listed. WP | WP name | Responsible institution | WP Reference
---|---|---|---
1 | PSF-R management and systems engineering | INAF-OAPD |
2 | PSF-R SCAO algorithm | JKU Linz | [13]
3 | PSF-R MCAO algorithm | JKU Linz | [13]
4 | SCAO simulations for PSF-R | LESIA | [12, 14]
5 | MCAO simulations for PSF-R | INAF-OAPD | [15, 16]
6 | Non-AO PSF-R algorithms and calibrations | INAF-OAR |
7 | PSF-R data flow architecture | NOVA | [17, 18]
8 | PSF-R user software | INAF-OAR |
9 | PSF-R scientific evaluation | INAF-OAPD | [19, 20]

The WP1 (PSF-R management and systems engineering) is responsible for the
communication within the PSF-R team and coordinating the project. It is in
charge of the managerial and sociological aspects within the team, of the
development plan and its coordination with the rest of the MICADO and MAORY
schedule.

The algorithms for the reconstruction of the SCAO PSF are developed within WP2
(PSF-R SCAO algorithm). The data items, that need to be provided by the Real-
Time Control (RTC), and the relevant calibrations are defined, in coordination
with the MICADO SCAO team. Test cases provided by the SCAO simulation WP (WP4)
are used to check the performance of the SCAO algorithms. Splitting
simulations and reconstruction tasks avoids development biases. The error
budget is estimated by measuring the noise of the input data and propagating
the errors throughout the algorithms. This WP contributes the SCAO
reconstruction module to the final PSF-R service.

Similarly, WP3 (PSF-R MCAO algorithm) provides the algorithms to reconstruct
the PSF for the MCAO configuration with MAORY. The development will be done in
close collaboration with the MAORY team. It should be noticed that MAORY
follows a different schedule than MICADO, which might have an impact on the
work planning.

As already mentioned, the simulations of the AO systems should follow as close
as possible the development of the systems themselves. Therefore, WP4 (SCAO
simulations for PSF-R) is split from the reconstruction algorithms and
assigned to the team building the AO system. In this way the in-depth
knowledge of the instrument development team can be used and increases the
confidence that the PSF reconstruction will perform as expected on the final
system. Additionally, the development of the SCAO system is supported by a
dedicated optical bench breadboard and an on-sky test. Both these experiments
can improve the insight on the PSF reconstruction performance and the
technical readiness of the algorithms. So it is crucial for the PSF-R
development to be part of those tests. The WP4 is waiting for test requests
from the SCAO algorithm development team.

For WP5 (MCAO simulations for PSF-R) similar considerations apply as in the
case of WP4, so it will be assigned to the MAORY team. There is no dedicated
bench experiment foreseen in the MAORY development at the moment, but there is
the possibility for a test during the integration phase before shipping MAORY
to Chile (planned for 2025). Also WP5 is waiting for test requests from the
MCAO algorithm development team.

There are several systems influencing the PSF beside the classically
considered AO system. The WP6 (Non-AO PSF-R algorithms and calibration) is
concerned with these other impacts on the PSF and investigates their
influence, the data available and the calibrations necessary for compensation.
With an error budget all the terms are balanced and the work on the specific
elements prioritised. The output algorithms of this work package are either
used as input for the AO PSF algorithms or to post-process the PSF. The
necessary data items for the algorithms will be defined and communicated to
the WP7 for proper collection.

The WP7 (PSF-R data flow architecture) works on the design of the data flow
system that starts from the raw data in the various systems, as defined by the
algorithmic work packages, and ends with the storage in the PSF-R data
archive. This work and the subsequent qualification of the implementation is
performed in close collaboration with ESO.

The WP8 (PSF-R user software) will integrate the PSF-R modules and algorithms
developed in the algorithmic work packages in a user software. The software
has to be written and documented as agreed by ESO so that it can be accepted
for takeover.

Finally, the goal of the WP9 (PSF-R scientific evaluation) is to provide the
algorithmic development (the WP2 and WP3) with performance metrics, define
science test cases and evaluate the impact of the PSF reconstruction
algorithms on the science test cases.

##  3 PSF RECONSTRUCTION

###  3.1 Why the PSF-R is needed

Given the unprecedented size of the ELT, the resulting PSF is quantitatively
different from that of any present AO facility. In fact, the core of the ELT
PSF is a factor $\sim 5$ narrower in diameter. This means that with respect to
the core, and at the same Strehl ratio, (i) the surface brightness of the halo
is a factor $\sim 5^{2}$ fainter, and (ii) the size of the halo is a factor
$\sim 5$ larger. This means that, for most stars the halo is barely detected
as a distinct feature, and that multiple halos blend into a faint distributed
background. On the other hand, the first Airy ring will be a prominent feature
of the PSF halo.

These conditions make the measure of the total flux of single sources
particularly challenging, without external input (like e.g. an indication of
the Strehl ratio), consequently affecting photometric measurements. The lower
is the achieved Strehl ratio of the observation, the higher is the
contribution of the PSF halo to the total flux of the source.

The PSF-R philosophy adopted for MICADO is to reconstruct a PSF independently
of the science data, to have the widest applicability (see Wagner et al.
2019[21] for a review on PSF determination techniques). In fact, it has been
estimated, referring to extragalactic targets, that in the $35$-$40\%$ of
cases, a MICADO science frame is void of point sources suitable for standard
PSF characterization[22] . In MICADO PSF-R, the AO telemetry is used to model
the optical transfer function (OTF), together with nearby $C_{n}^{2}$ profiler
data, and on-sky phase diversity. Specifically, an initial PSF estimate is
adapted according to various degradation effects and other calibrations, to
provide a final reconstructed PSF. It is important to handle separately the AO
and non-AO effects, as their scaling with wavelength is different. The
residual wavefront leads to a wavelength dependence of the PSF shape; residual
chromatism affects the location of the PSF at each wavelength; and residual
vibration is independent of wavelength. The plan for MICADO is that these
effects are calibrated or measured so that they can be included appropriately
in the final reconstructed PSF. Crucially, the data flow architecture has to
be compatible with ESO’s data-flow and archiving systems, especially on how
the various data are gathered, saved, and made available to the user. In
addition, a significant effort is being made to derive non-AO contributions in
a way that minimizes on-sky calibrations.

###  3.2 Algorithm for PSF-R

The PSF obtained by the optical complex system made of the chain of the
turbulent atmosphere, the ELT telescope system, the MICADO relay optics, the
SCAO system and the MICADO science instrument will not be perfect. Along with
the instrument specifications, the SCAO system will deliver a PSF whose Strehl
ratio is up to $60\%$ (Goal: $70\%$) in K-band in median atmospheric
turbulence conditions with a bright NGS in the close vicinity of the science
target. The optical performance will be reduced in less optimal cases (strong
turbulence, less bright NGS, science object far away from the science target).

PSF-Reconstruction algorithms mostly operate to produce an OTF (i.e. the
Fourier space counterpart of the PSF). The OTF of the PSF part related to
atmospheric turbulence in a specific direction and wavelength is computed to
reconstruct the post-AO PSF associated with a given MICADO SCAO scientific
observation. The OTF can be calculated from the so-called structure function
($D$) which is related to the residual incoming phase. The approach from Veran
et al. (1997)[23] is followed, splitting the structure function into several
parts which are assumed to be independent. Each of these parts will be
wavelength dependent and use calibration data items. The final result is the
combined post-AO OTF for a specific wavelength $\lambda$ and direction
$\alpha$ associated to one observational frame. It needs to be multiplied by
the non-AO OTF and Fourier transformed to obtain the PSF. All computations are
performed on the same resolution otherwise interpolations will add ghosts when
applying a Fourier Transform. As regular grids for sampling the incoming phase
are used, particular care is given in making sure that with a hexagonal
actuator pattern, all frequencies which can be represented on M4 are still
sampled. Using the currently available ESO data package, this need is
fulfilled with a grid spacing of $0.25$ m projected on M1. Furthermore, the
phase $\phi$ is related to the wavefront $\varphi$ via
$\phi=\frac{2\pi}{\lambda}\varphi$ and thus wavelength dependent.

The algorithm developed in the context of the MICADO PSF-R is currently
composed of the following blocks:

  1. 1.

Wavefront low order structure function. From the reconstructed residual
incoming phase $\phi_{rec,t}$, for time frames $t$ and related to the
direction of the NGS at wavelength $\lambda$, the 2D averaged structure
function $D_{rec}$ is computed. $\phi_{rec,t}$ is computed from the WFS
measurements using the AO control algorithm. Note that this computation can be
done as a pre-processing step. Alternatively the computed deformable mirrors
(DM) shapes (i.e. M4 + M5 commands in SCAO) can be used to extract this
information.

  2. 2.

Aliasing structure function. A model is used to propagate simulated higher
order parts of phases through the WFS and then compute the response of the
control algorithm. The output is the structure function of the aliasing part
of the phase $D_{al}$. This procedure is rerun for multiple time steps before
all these computed aliasing phases are used in a structure function
computation. The input list for this algorithm includes seeing, $r_{0}$, used
AO control matrix/algorithm and Pupil mask $P$.

  3. 3.

Noise structure function. The AO control algorithm is used to propagate the
noise model from the WFS level to wavefront level $D_{n}$. This computation
needs the following inputs: guide star magnitude (or photon flux level) at
$\lambda$; Pupil mask $P$ (from calibrations); M4 actuator positions;
deformable mirrors influence functions and active WFS subapertures.

  4. 4.

Wavefront higher order structure function. From a statistical model of the
atmospheric turbulence (von Karman, Kolmogorov) a structure function for the
higher order part not corrected by the AO system, $D_{\perp}$, is computed.
This cannot be done once for all since the pupil mask changes.

  5. 5.

Anisoplanatic structure function. Using the reconstructed residual incoming
phases $\phi_{t,rec}$, and knowledge on the atmospheric turbulence profile, an
atmospheric tomography can be performed to recover the atmospheric turbulence
residual in the direction of interest. The difference between this residual
and the reconstructed residual wavefronts during the AO run, the anisoplanatic
structure function $D_{an}$ can be computed.

  6. 6.

Combined post-AO OTF. From each of the above structure functions ($D_{rec}$,
$D_{al}$, $D_{n}$, $D_{\perp}$ and $D_{an}$) an OTF has to be computed as

$OTF_{x}=\exp(-\frac{1}{2}D_{x})$

with $x\in\\{\|,\perp,an\\}$, where $D_{\|}=D_{rec}-D_{n}+D_{al}$. The
telescope OTF ($OTF_{tel}$), instead, is computed from $P$ and these OTFs are
combined by multiplication to obtain the full AO OTF:

$OTF_{AO}=OTF_{tel}\cdot OTF_{\|}\cdot OTF_{\perp}\cdot OTF_{an}$

.

Finally, the Non-AO effects are also added. They are composed by two parts:
the non-common path aberrations (NCPA) and the vibration pattern. The NCPA are
a series of Zernike coefficients ($Z_{3},...,Z_{100}$) that are transformed
into an OTF function defined as $OTF_{NCPA}$: they model the NCPA for that
specific instrument configuration and telescope pointing. The vibration
pattern is a 2D array where each element is the integral over the exposure
time of the presence of the chief ray in that element, and the whole array is
normalised to $total(array)=1$: in other words, it is equivalent to a frame in
which the value of each pixel is proportional to the number of times the chief
ray has passed over it. The final PSF, $PSF_{MICADO}$, is retrieved in three
steps for each point of the field:

  1. 1.

Compute $OTF_{NCPA}$ from Zernike coefficients.

  2. 2.

The $OTF_{AO}$ is combined with the $OTF_{NCPA}$ to calculate the
$OTF_{AO+NCPA}$ as

$OTF_{AO+NCPA}=OTF_{AO}\cdot OTF_{NCPA}$

  3. 3.

The $OTF_{AO+NCPA}$ is transformed via $FFT^{-1}$ into the $PSF_{AO+NCPA}$

  4. 4.

The kernel of the residual vibrations pattern $K_{vibration}$ is registered to
the focal plane coordinates, and convolved with the $PSF_{AO+NCPA}$ to obtain
the final $PSF_{MICADO}$ as

$PSF_{MICADO}=PSF_{AO+NCPA}*K_{vibration}$

##  4 EXAMPLE OF PSF RECONSTRUCTION

The first prototype of the PSF-R algorithm has been tested using closed loop
telemetry and near infrared images acquired with the SOUL[24] system mounted
on the front bent Gregorian focus of the LBT[25] telescope and the LUCI[26]
camera. A diffraction limited camera (a.k.a. N30) equips the LUCI instrument
for both imaging and spectroscopy, offering a $2k\times 2k$ sensitive area
corresponding to $30"\times 30"$ on the sky (14.95 mas/px). For this
preliminary testing of the PSF-R algorithm, laboratory-like condition has been
used: optical turbulence disturbance is generated and corrected[27] using the
Adaptive Secondary AdSec[28] of the LBT, with the telescope and AO system
configured for daytime[29]. The SOUL team provided the essential data set for
the PSF-R: AO WFS slopes data history, control matrix, interaction matrix,
gain vector, and the pupil definition. The PSF-R does not use the camera
observing data for the analysis. The LUCI observed PSF was only used as
comparison for the analysis of the performed PSF-R. Figure 1 shows the
matching of the observed and reconstructed PSF. The central bright area
corresponds to the spatial frequencies controlled by the AdSec and limited by
inter-actuator distance. Beyond, the uncorrected dome seeing turbulence
contributes are negligible. More quantitatively, in Figure 2 it can be noticed
the very good correspondence of the measured and reconstructed Strehl ratio
(FeII filter) of 0.58 and 0.61 respectively. The PSF-R algorithm also
reconstructed the FWHM and encircled energy profiles within the measurements
errors (within a few percent points).

Figure 1: The picture shows: (Left) the LUCI observed PSF, (Right) the
reconstructed PSF. The frame angular size is $3".8\times 3".8$.

Figure 2: Intensity and encircled energy profiles.

##  5 SUMMARY

The PSF-R capability is a deliverable of the MICADO@ESO-ELT project. This will
make MICADO the first instrument to date to offer this service. In order for
the deliverable to have the widest scientific applicability, a pure PSF-R
approach is pursued, i.e. to reconstruct the PSF independently of the science
data. The MICADO PSF-R team is working on the development of the PSF-R method
with $9$ WPs actively interacting. Their aims and responsibilites have been
presented, along with a description of the PSF-R algorithm. The main features
of the MICADO PSF-R service include:

  * –

Generate field-dependent reconstructed PSF at the desired wavelength, taking
into account both AO and non-AO (e.g. telescope or atmospheric) input.

  * –

Provide reconstructed PSFs through an archive querying system to restore the
telemetry data synchronous to each of the science frames that MICADO will
generate.

  * –

Produce the output according to user specifications.

  * –

Support the state-of-the-art scientific analysis of the MICADO imaging and
spectroscopic data (both for SCAO and MCAO observing mode, the latter provided
by the MAORY module).

Finally, results coming from the application of the first prototype of the
PSF-R algorithms on real SOUL@LBT data have been presented. The obtained
reconstructed PSF shows small deviation from the observed one, of the order of
few percent points in Strehl ratio, FWHM and encircle energy profile.
Remaining within the measurement errors, these results are encouraging already
at the present stage and they are a good starting point toward a successful
development of the PSF-R algorithm for MICADO. The PSF-R team of MICADO will
work in the next $\sim 5$ years to further improve these algorithms. The
general philosophy of MICADO PSF-R could be exported also to other AO
facilities of ESO-ELT and other gigantic telescopes in the future.

###### Acknowledgements.

This work has been partly supported by INAF through the Math, ASTronomy and
Research (MAST&R), a working group for mathematical methods for high-
resolution imaging. The PSF-R Team thanks Miska Le Louarn (ESO) for continued
support and enthusiasm for the project.

## References

  * [1] Davies, R., Ageorges, N., Barl, L., Bedin, L. R., Bender, R., Bernardi, P., Chapron, F., Clenet, Y., Deep, A., Deul, E., Drost, M., Eisenhauer, F., Falomo, R., Fiorentino, G., Förster Schreiber, N. M., Gendron, E., Genzel, R., Gratadour, D., Greggio, L., Grupp, F., Held, E., Herbst, T., Hess, H. J., Hubert, Z., Jahnke, K., Kuijken, K., Lutz, D., Magrin, D., Muschielok, B., Navarro, R., Noyola, E., Paumard, T., Piotto, G., Ragazzoni, R., Renzini, A., Rousset, G., Rix, H. W., Saglia, R., Tacconi, L., Thiel, M., Tolstoy, E., Trippe, S., Tromp, N., Valentijn, E. A., Verdoes Kleijn, G., and Wegner, M., “MICADO: the E-ELT adaptive optics imaging camera,” in [Proc. SPIE ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 7735, 77352A (Jul 2010).
  * [2] Gilmozzi, R. and Spyromilio, J., “The European Extremely Large Telescope (E-ELT),” The Messenger 127, 11 (2007).
  * [3] Gilmozzi, R. and Spyromilio, J., “The 42m European ELT: status,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 7012 (2008).
  * [4] Gilmozzi, R. and Spyromilio, J., “The 42m European ELT: status,” in [Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series ], Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 7012 (2008).
  * [5] Babcock, H. W., “Adaptive optics revisited,” Science 249(4966), 253–257 (1990).
  * [6] Beckers, J. M., “Increasing the size of the isoplanatic patch with multiconjugate adaptive optics.,” in [ESO Conference on Very Large Telescopes and their Instrumentation ], 2, 693–703 (1988).
  * [7] Beckers, J. M., “Detailed compensation of atmospheric seeing using multiconjugate adaptive optics.,” in [Active Telescope Systems ], Proc. SPIE 1114, 215–217 (1989).
  * [8] Ragazzoni, R., “Adaptive optics for giant telescopes: NGS vs. LGS,” in [Proceedings of the Backaskog workshop on extremely large telescopes ], Andersen, T., Ardeberg, A., and Gilmozzi, R., eds., 175–180 (2000).
  * [9] Ragazzoni, R., Marchetti, E., and Valente, G., “Adaptive-optics corrections available for the whole sky,” Nature 403, 54–56 (Jan. 2000).
  * [10] Ragazzoni, R., Marchetti, E., and Rigaut, F., “Modal tomography for adaptive optics,” A&A 342, L53–L56 (Feb. 1999).
  * [11] Diolaiti, E., “MAORY: A Multi-conjugate Adaptive Optics RelaY for the E-ELT,” The Messenger 140, 28–29 (2010).
  * [12] Clenet, Y., Esposito, S., Buey, T., Riccardi, A., Rousset, G., Lombini, M., Cohen, M., Spano, P., Gendron, E., Diolaiti, E., and Davies, R., “Implementation of SCAO for ELT-CAM / MICADO-MAORY,” (2015).  Accepted: 2020-03-20T10:48:56Z.
  * [13] Wagner, R., Hofer, C., and Ramlau, R., “Point spread function reconstruction for single-conjugate adaptive optics on extremely large telescopes,” Journal of Astronomical Telescopes, Instruments, and Systems 4, 049003 (Oct. 2018).
  * [14] Clénet, Y., Buey, T., Gendron, E., Hubert, Z., Vidal, F., Cohen, M., Chapron, F., Sevin, A., Fédou, P., Barbary, G., Baudoz, P., Borgo, B., Ben Nejma, S., Chambouleyron, V., Déo, V., Dupuis, O., Durand, S., Ferreira, F., Gaudemard, J., Gratadour, D., Huby, E., Huet, J.-M., Le Ruyet, B., Nguyen-Tuong, N., Perrot, C., Thijs, S., Younès, Y., Rousset, G., Feautrier, P., Zins, G., Diolaiti, E., Ciliegi, P., Esposito, S., Busoni, L., Schubert, J., Hartl, M., Hörmann, V., and Davies, R., “The MICADO first-light imager for the ELT: towards the preliminary design review of the MICADO-MAORY SCAO,” in [Adaptive Optics Systems VI ], Close, L. M., Schreiber, L., and Schmidt, D., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 1070313 (July 2018).
  * [15] Arcidiacono, C., Schreiber, L., Bregoli, G., Foppiani, I., Oberti, S., Diolaiti, E., Agapito, G., Puglisi, A., Xompero, M., Cortecchia, F., Patti, M., Lombini, M., Busoni, L., Vérinaud, C., Felini, C., De Caprio, V., Cosentino, G., Ciliegi, P., Bellazzini, M., Feautrier, P., Esposito, S., and Ragazzoni, R., “Numerical simulations of MAORY MCAO module for the ELT,” in [Adaptive Optics Systems VI ], Close, L. M., Schreiber, L., and Schmidt, D., eds., Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series 10703, 107034I (July 2018).
  * [16] Portaluri, E., Viotto, V., Ragazzoni, R., Arcidiacono, C., Bergomi, M., Dima, M., Greggio, D., Farinato, J., and Magrin, D., “Simulations of ELT-GMCAO Performance for Deep Field Observations,” PASP 132, 084502 (Aug. 2020).
  * [17] Verdoes Kleijn, G., “Science Data Management for the E-ELT: usecase MICADO,” in [Science Operations 2015: Science Data Management ], 1 (Dec. 2015).
  * [18] Buddelmeijer, H., Verdoes Kleijn, G. A., and Leschinski, K., “Data model as agile basis for evolving calibration software,” arXiv e-prints arXiv:2011.04516, arXiv:2011.04516 (Nov. 2020).
  * [19] Gullieuszik, M., Greggio, L., Falomo, R., Schreiber, L., and Uslenghi, M., “Probing the nuclear star cluster of galaxies with extremely large telescopes,” A&A 568, A89 (Aug. 2014).
  * [20] Gullieuszik, M., Falomo, R., Greggio, L., Uslenghi, M., and Fantinel, D., “Exploring high-z galaxies with the E-ELT,” A&A 593, A24 (Aug. 2016).
  * [21] Wagner, R., Beltramo-Martin, O., Correia, C., Fétick, R., Ramlau, R., Fusco, T., and Neichel, B., “Overview of PSF determination techniques for adaptive-optics assisted ELT instruments,” in [6th International Conference on Adaptive Optics for Extremely Large Telescopes, AO4ELT 2019 ], (June 2019).
  * [22] Davies, R., “An Astrophysical Context for PSF Reconstruction,” tech. rep., MPE - Max-Planck-Institut für extraterrestrische Physik (09 2019).
  * [23] Veran, J. P., Rigaut, F., Maitre, H., and Rouan, D., “Estimation of the adaptive optics long-exposure point-spread function using control loop data.,” Journal of the Optical Society of America A 14, 3057–3069 (Nov. 1997).
  * [24] Pinna, E., Esposito, S., Hinz, P., Agapito, G., Bonaglia, M., Puglisi, A., Xompero, M., Riccardi, A., Briguglio, R., Arcidiacono, C., Carbonaro, L., Fini, L., Montoya, M., and Durney, O., “SOUL: the Single conjugated adaptive Optics Upgrade for LBT,” in [Adaptive Optics Systems V ], 9909, 99093V, International Society for Optics and Photonics (July 2016).
  * [25] Hill, J. M., “The Large Binocular Telescope,” Applied Optics 49, 115–122 (Apr. 2010).
  * [26] Seifert, W., Appenzeller, I., Baumeister, H., Bizenberger, P., Bomans, D., Dettmar, R.-J., Grimm, B., Herbst, T., Hofmann, R., Juette, M., Laun, W., Lehmitz, M., Lemke, R., Lenzen, R., Mandel, H., Polsterer, K., Rohloff, R.-R., Schuetze, A., Seltmann, A., Thatte, N. A., Weiser, P., and Xu, W., “LUCIFER: a multimode NIR instrument for the LBT,” in [Instrument Design and Performance for Optical/Infrared Ground-based Telescopes ], 4841, 962–973, International Society for Optics and Photonics (Mar. 2003).
  * [27] Esposito, S., Riccardi, R., Fini, L., Puglisi, A. T., Pinna, E., Xompero, M., Briguglio, R., Quirós-Pacheco, F., Stefanini, P., Guerra, J. C., Busoni, L., Tozzi, A., Pieralli, F., Agapito, G., Brusa-Zappellini, G., Demers, R., Brynnel, J., Arcidiacono, C., and Salinari, P., “First light AO (FLAO) system for LBT: final integration, acceptance test in Europe, and preliminary on-sky commissioning results,” in [Adaptive Optics Systems II ], Proc. SPIE 7736 (2010).
  * [28] Riccardi, A., Brusa, G., Salinari, P., Gallieni, D., Biasi, R., Andrighettoni, M., and Martin, H. M., “Adaptive secondary mirrors for the Large Binocular Telescope,” Proceedings of SPIE 4839, 721–732 (2003).
  * [29] Riccardi, A., Xompero, M., Zanotti, D., Busoni, L., Vecchio, C. D., Salinari, P., Ranfagni, P., Zappellini, G. B., Biasi, R., Andrighettoni, M., Gallieni, D., Anaclerio, E., Martin, H. M., and Miller, S. M., “Adaptive secondary mirror for the Large Binocular Telescope: results of acceptance laboratory test,” Proceedings of SPIE 7015(37), 701512 (2008).


# Direct constructions of some group divisible designs with block size $4$ and
up to $50$ points

R. Julian R. Abel
School of Mathematics and Statistics
UNSW Sydney
NSW 2052, Australia
<EMAIL_ADDRESS>Thomas Britz
School of Mathematics and Statistics
UNSW Sydney
NSW 2052, Australia
<EMAIL_ADDRESS>Yudhistira A. Bunjamin
School of Mathematics and Statistics
UNSW Sydney
NSW 2052, Australia
<EMAIL_ADDRESS>Diana Combe
School of Mathematics and Statistics
UNSW Sydney
NSW 2052, Australia
<EMAIL_ADDRESS>
Abstract: In this note, we give direct constructions of some group divisible
designs (GDDs) with block size $4$ that have up to $50$ points.

Keywords: group divisible design (GDD), feasible group type.

Mathematics Subject Classification: 05B05

##  1 Introduction

Let $X$ be a finite set of points with a partition into parts which we call
groups. Any $k$-element subset of $X$ is called a block. A collection of
blocks is a group divisible design with block size $k$, or $k$-GDD if (i) no
two points from the same group appear together in any block and (ii) any two
points from distinct groups appear together in exactly one block. The group
type (or type) of a $k$-GDD is the multiset $\left\\{|G|:G{\text{ is a part of
the partition}}\right\\}$. The group type can also be expressed in
‘exponential’ notation where the type $t_{1}^{u_{1}}t_{2}^{u_{2}}\ldots
t_{m}^{u_{m}}$ means there are $u_{i}$ groups of size $t_{i}$ for
$i=1,2,\dots,m$.

There are known necessary conditions for the existence of a $4$-GDD of type
$\\{g_{1},g_{2},\ldots,g_{m}\\}$. These are given in [4]. However, these
necessary conditions are not sufficient.

In this note, we are concerned with $4$-GDDs. We say that a multiset
$\\{g_{1},g_{2},\ldots,g_{m}\\}$ of positive integers is a feasible group type
for a $4$-GDD if it satisfies the known necessary conditions.

In [5], Kreher and Stinson gave the list of feasible group types for $4$-GDDs
with up to $v=30$ points. There the existence question was answered for all
but three group types. The existence results for those three types were
completed in [1, 3].

In [4], the list of feasible group types for $4$-GDDs was extended to those
with up to $v=50$ points. For $31\leq v\leq 50$ and $v\equiv 0\pmod{3}$ there
exist $4$-GDDs for all feasible group types; this result was completed in [4].
In [2], the feasible types for $31\leq v\leq 50$ and $v\equiv 1$ and
$2\pmod{3}$ were considered. The questions of existence were completed for
$v\equiv 1\pmod{3}$; and for $v\equiv 2\pmod{3}$ the known results were
extended leaving unknown the question of existence for types
$2^{11}8^{1}11^{1}$, $2^{1}5^{4}8^{1}11^{1}$, $2^{6}5^{2}11^{2}$,
$2^{5}5^{3}8^{1}11^{1}$, $2^{2}5^{2}8^{1}11^{2}$, $2^{1}5^{3}8^{2}11^{1}$,
$2^{5}5^{3}11^{2}$, $2^{2}5^{2}11^{3}$, $2^{1}5^{3}8^{1}11^{2}$,
$5^{4}8^{2}11^{1}$, $2^{9}5^{2}11^{2}$, $2^{8}5^{3}8^{1}11^{1}$,
$2^{6}5^{1}11^{3}$, $2^{5}5^{2}8^{1}11^{2}$, $2^{4}5^{4}8^{1}14^{1}$,
$2^{4}5^{3}8^{2}11^{1}$, $2^{3}11^{4}$, $2^{2}5^{7}11^{1}$,
$2^{2}5^{1}8^{1}11^{3}$, $2^{1}5^{2}8^{2}11^{2}$ and $5^{3}8^{3}11^{1}$.

In this note, we provide direct constructions of $4$-GDDs with up to $50$
points, particularly for the group types for which the question of existence
was left open in [2]. The primary intention of this note is to make these
constructions publicly available before they appear in future publications by
the current authors. Occasionally, we may also include constructions of
$4$-GDDs of group types for which existence has already been determined.
However, these constructions will be different from the ones that were first
used to prove the existence of $4$-GDDs of these group types. As a result,
this note may be updated from time to time with new constructions.

##  2 Direct constructions by assuming an automorphism

The $4$-GDDs in this section have been found directly by assuming the
existence of a cyclic automorphism group, $\mathbb{G}$. For each $4$-GDD, a
set of base blocks is given. The blocks are obtained by developing the
subscripts of the points in each base block from copies of $\mathbb{G}$ over
$\mathbb{G}$; the infinite points remain unaltered when developed.

We note that a $4$-GDD of type $3^{8}6^{1}12^{1}$ was already known to exist
prior to this note. The $4$-GDD of type $3^{8}6^{1}12^{1}$ in [6] was obtained
by Rees and Stinson implicitly. A $\\{3,4\\}$-GDD is defined similarly to a
$k$-GDD except that each block is either of size $3$ or size $4$. Rees and
Stinson construct a $\\{3,4\\}$-GDD of type $3^{8}6^{1}$ whose triples fall
into $12$ resolution classes in [6, Appendix]. To construct a $4$-GDD of type
$3^{8}6^{1}12^{1}$, they started with the $\\{3,4\\}$-GDD of type $3^{8}6^{1}$
and turned the blocks of size $3$ into blocks of size $4$ by adding each point
from the group of size $12$ to the blocks of size $3$ in one of the resolution
classes.

In Example 2.1, we give a $4$-GDD of type $3^{8}6^{1}12^{1}$ that was
constructed directly.

######  Example 2.1.

A $4$-GDD of type $3^{8}6^{1}12^{1}$ can be constructed in the following
manner.

The points in the $4$-GDD are

  * •

$a_{i},b_{i},q_{i},r_{i},y_{i},z_{i}$ for $i\in{\mathbb{Z}}_{6};$

  * •

$p_{i}$ for $i\in{\mathbb{Z}}_{3};$ and

  * •

$\infty_{1},\infty_{2},\infty_{3}$.

The groups are

  * •

$\\{a_{i}:i\in{\mathbb{Z}}_{6}\\}\cup\\{b_{i}:i\in{\mathbb{Z}}_{6}\\};$

  * •

$\\{p_{i}:i\in{\mathbb{Z}}_{3}\\}\cup\\{\infty_{1},\infty_{2},\infty_{3}\\};$

  * •

$\\{q_{i},q_{i+2},q_{i+4}\\}$, $\\{r_{i},r_{i+2},r_{i+4}\\}$,
$\\{y_{i},y_{i+2},y_{i+4}\\};$ and

  * •

$\\{z_{i},z_{i+2},z_{i+4}\\}$ for $i\in\\{0,1\\}$.

Blocks are obtained by developing modulo $6$ the subscripts of the noninfinite
points in the base blocks given in Table 1. The first two blocks in the first
column generate three blocks each.

Table 1: Base blocks for the $4$-GDD of type $3^{8}6^{1}12^{1}$ in Example 2.1
$\begin{array}[]{|l|l|l|l|l|l|}\hline\cr\\{p_{0},p_{3},q_{0},q_{3}\\}&\\{a_{0},p_{0},p_{1},y_{0}\\}&\\{a_{0},q_{3},u_{5},\infty_{3}\\}&\\{b_{0},p_{4},q_{5},u_{0}\\}&\\{b_{0},p_{5},u_{4},u_{5}\\}\\\
\hline\cr\\{t_{0},t_{3},u_{0},u_{3}\\}&\\{a_{0},p_{2},q_{0},t_{0}\\}&\\{a_{0},q_{1},t_{3},u_{1}\\}&\\{b_{0},p_{0},q_{2},\infty_{2}\\}&\\{b_{0},q_{0},q_{1},t_{4}\\}\\\
\hline\cr&\\{a_{0},p_{3},q_{2},u_{0}\\}&\\{a_{0},q_{4},t_{5},y_{2}\\}&\\{b_{0},p_{2},u_{3},y_{0}\\}&\\{b_{0},q_{4},t_{3},y_{1}\\}\\\
\hline\cr&\\{a_{0},p_{4},t_{1},\infty_{1}\\}&\\{a_{0},q_{5},u_{2},y_{1}\\}&\\{b_{0},p_{1},t_{1},t_{2}\\}&\\{b_{0},q_{3},u_{2},\infty_{1}\\}\\\
\hline\cr&\\{a_{0},p_{5},t_{4},u_{3}\\}&\\{a_{0},t_{2},u_{4},\infty_{2}\\}&\\{b_{0},p_{3},t_{5},\infty_{3}\\}&\\{b_{0},t_{0},u_{1},y_{2}\\}\\\
\hline\cr\end{array}$

######  Lemma 2.2.

There exists a $4$-GDD of type $2^{6}5^{2}11^{2}$.

###### Proof.

The points in the $4$-GDD are

  * •

$a_{i},b_{i},c_{i},f_{i},g_{i},p_{i},q_{i},r_{i},s_{i},t_{i}$ for
$i\in{\mathbb{Z}}_{4}$; and

  * •

$h_{i},u_{i}$ for $i\in{\mathbb{Z}}_{2}$.

The groups are

  * •

$\\{a_{i},b_{i}\\}$ for $i\in{\mathbb{Z}}_{4}$;

  * •

$\\{c_{i},c_{i+2}\\}$;

  * •

$\\{f_{i},f_{i+2},g_{i},g_{i+2},h_{i}\\}$; and

  * •

$\\{p_{i},p_{i+2},q_{i},q_{i+2},r_{i},r_{i+2},s_{i},s_{i+2},t_{i},t_{i+2},u_{i}\\}$
for $i\in\\{0,1\\}$.

Blocks are obtained by developing modulo $4$ the subscripts of the noninfinite
points in the base blocks given in Table 2. The first block in the first
column generates one block, while the second block in that column generates
two blocks. ∎

Table 2: Base blocks for the $4$-GDD of type $2^{6}5^{2}11^{2}$ in Lemma 2.2
$\begin{array}[]{|l|l|l|l|l|l|}\hline\cr\\{h_{0},h_{1},u_{0},u_{1}\\}&\\{a_{0},a_{1},p_{0},r_{1}\\}&\\{a_{0},g_{2},r_{3},s_{2}\\}&\\{b_{0},f_{3},t_{0},u_{1}\\}&\\{c_{0},g_{1},r_{3},u_{0}\\}\\\
\hline\cr\\{a_{0},a_{2},b_{1},b_{3}\\}&\\{a_{0},b_{2},c_{0},h_{0}\\}&\\{a_{0},g_{3},r_{2},t_{1}\\}&\\{b_{0},g_{2},p_{0},s_{3}\\}&\\{c_{0},h_{1},r_{2},s_{3}\\}\\\
\hline\cr&\\{a_{0},c_{1},f_{1},g_{0}\\}&\\{a_{0},h_{1},q_{1},t_{0}\\}&\\{b_{0},g_{3},p_{3},t_{2}\\}&\\{f_{0},f_{1},p_{3},s_{0}\\}\\\
\hline\cr&\\{a_{0},c_{2},p_{1},p_{2}\\}&\\{b_{0},b_{1},q_{2},s_{1}\\}&\\{b_{0},g_{0},r_{0},t_{1}\\}&\\{f_{0},h_{1},p_{0},r_{3}\\}\\\
\hline\cr&\\{a_{0},c_{3},s_{0},s_{3}\\}&\\{b_{0},c_{1},f_{0},g_{1}\\}&\\{b_{0},h_{1},s_{2},t_{3}\\}&\\{g_{0},g_{1},q_{2},s_{3}\\}\\\
\hline\cr&\\{a_{0},f_{2},q_{2},q_{3}\\}&\\{b_{0},c_{3},p_{1},u_{0}\\}&\\{c_{0},c_{1},q_{3},t_{0}\\}&\\{g_{0},h_{1},p_{1},q_{0}\\}\\\
\hline\cr&\\{a_{0},f_{3},s_{1},u_{0}\\}&\\{b_{0},c_{0},q_{0},r_{1}\\}&\\{c_{0},f_{2},q_{1},r_{0}\\}&\\\
\hline\cr&\\{a_{0},f_{0},t_{2},t_{3}\\}&\\{b_{0},f_{1},p_{2},q_{3}\\}&\\{c_{0},f_{1},s_{2},t_{1}\\}&\\\
\hline\cr&\\{a_{0},g_{1},q_{0},u_{1}\\}&\\{b_{0},f_{2},r_{2},r_{3}\\}&\\{c_{0},g_{2},p_{1},t_{2}\\}&\\\
\hline\cr\end{array}$

######  Lemma 2.3.

There exists a $4$-GDD of type $2^{9}5^{2}11^{2}$.

###### Proof.

The points in the $4$-GDD are

  * •

$a_{i},b_{i},c_{i},d_{i},f_{i},g_{i},p_{i},q_{i},r_{i},s_{i},t_{i}$ for
$i\in{\mathbb{Z}}_{4}$;

  * •

$h_{i},u_{i}$ for $i\in{\mathbb{Z}}_{2}$; and

  * •

$\infty_{1},\infty_{2}$.

The groups are

  * •

$\\{a_{i},b_{i}\\}$ and $\\{c_{i},d_{i}\\}$ for $i\in{\mathbb{Z}}_{4}$;

  * •

$\\{f_{i},f_{i+2},g_{i},g_{i+2},h_{i}\\}$;

  * •

$\\{p_{i},p_{i+2},q_{i},q_{i+2},r_{i},r_{i+2},s_{i},s_{i+2},t_{i},t_{i+2},u_{i}\\}$
for $i\in\\{0,1\\}$; and

  * •

$\\{\infty_{1},\infty_{2}\\}$.

Blocks are obtained by developing modulo $4$ the subscripts of the noninfinite
points in the base blocks given in Table 3. The first block in the first
column generates one block, while the other four blocks in that column
generate two blocks each. ∎

Table 3: Base blocks for the $4$-GDD of type $2^{9}5^{2}11^{2}$ in Lemma 2.3
$\begin{array}[]{|l|l|l|l|l|}\hline\cr\\{h_{0},h_{1},u_{0},u_{1}\\}&\\{a_{0},a_{1},c_{0},g_{0}\\}&\\{a_{0},g_{1},p_{2},r_{1}\\}&\\{b_{0},g_{2},g_{3},q_{2}\\}&\\{d_{0},f_{1},p_{0},u_{1}\\}\\\
\hline\cr\\{a_{0},a_{2},u_{0},\infty_{1}\\}&\\{a_{0},b_{1},f_{0},h_{1}\\}&\\{a_{0},h_{0},q_{2},s_{3}\\}&\\{b_{0},g_{1},t_{1},u_{0}\\}&\\{d_{0},f_{2},r_{1},\infty_{2}\\}\\\
\hline\cr\\{b_{0},b_{2},u_{1},\infty_{2}\\}&\\{a_{0},b_{2},p_{1},s_{0}\\}&\\{a_{0},q_{3},s_{2},\infty_{2}\\}&\\{b_{0},h_{1},p_{2},r_{3}\\}&\\{d_{0},g_{2},q_{3},u_{0}\\}\\\
\hline\cr\\{c_{0},c_{2},h_{0},\infty_{2}\\}&\\{a_{0},b_{3},q_{0},t_{1}\\}&\\{b_{0},b_{1},c_{0},p_{1}\\}&\\{c_{0},c_{1},f_{0},q_{0}\\}&\\{d_{0},g_{0},s_{1},s_{2}\\}\\\
\hline\cr\\{d_{0},d_{2},h_{1},\infty_{1}\\}&\\{a_{0},c_{1},r_{0},s_{1}\\}&\\{b_{0},c_{1},d_{0},t_{0}\\}&\\{c_{0},d_{2},g_{1},r_{2}\\}&\\{d_{0},h_{0},p_{2},q_{1}\\}\\\
\hline\cr&\\{a_{0},c_{2},r_{2},u_{1}\\}&\\{b_{0},c_{2},g_{0},s_{0}\\}&\\{c_{0},d_{1},q_{1},t_{0}\\}&\\{f_{0},f_{1},s_{2},t_{3}\\}\\\
\hline\cr&\\{a_{0},d_{1},f_{1},g_{2}\\}&\\{b_{0},d_{1},d_{2},s_{1}\\}&\\{c_{0},f_{2},p_{0},p_{3}\\}&\\{f_{0},g_{3},q_{1},r_{2}\\}\\\
\hline\cr&\\{a_{0},d_{0},f_{3},p_{3}\\}&\\{b_{0},d_{3},r_{1},r_{2}\\}&\\{c_{0},f_{1},s_{1},u_{0}\\}&\\{g_{0},h_{1},r_{2},t_{3}\\}\\\
\hline\cr&\\{a_{0},d_{3},p_{0},q_{1}\\}&\\{b_{0},f_{1},q_{0},q_{3}\\}&\\{c_{0},g_{3},p_{2},t_{1}&\\{g_{0},p_{2},s_{3},\infty_{1}\\}\\\
\hline\cr&\\{a_{0},d_{2},t_{0},t_{3}\\}&\\{b_{0},f_{0},r_{0},s_{3}\\}&\\{c_{0},h_{1},s_{3},t_{2}\\}&\\{g_{0},p_{0},t_{1},\infty_{2}\\}\\\
\hline\cr&\\{a_{0},f_{2},r_{3},t_{2}\\}&\\{b_{0},f_{2},t_{3},\infty_{1}\\}&\\{c_{0},q_{2},r_{1},\infty_{1}\\}&\\\
\hline\cr\end{array}$

######  Lemma 2.4.

There exists a $4$-GDD of type $2^{6}5^{1}11^{3}$.

###### Proof.

The points in the $4$-GDD are

  * •

$a_{i},b_{i},c_{i},d_{i},e_{i},s_{i},t_{i}$ for $i\in{\mathbb{Z}}_{6}$;

  * •

$f_{i}$ for $i\in{\mathbb{Z}}_{3}$;

  * •

$p_{i},q_{i}$ for $i\in{\mathbb{Z}}_{2}$; and

  * •

$\infty$.

The groups are

  * •

$\\{a_{i},a_{i+3},b_{i},b_{i+3},c_{i},c_{i+3},d_{i},d_{i+3},e_{i},e_{i+3},f_{i}\\}$
for $i\in\\{0,1,2\\}$;

  * •

$\\{s_{i},t_{i}\\}$ for $i\in\mathbb{Z}_{6}$; and

  * •

$\\{p_{i},q_{i}\>:\>i\in{\mathbb{Z}}_{2}\\}\cup\\{\infty\\}$.

Blocks are obtained by developing modulo $6$ the subscripts of the noninfinite
points in the base blocks given in Table 4. The first three blocks in the
first column generate two blocks each, while the last two blocks in that
column generate three blocks each. ∎

Table 4: Base blocks for the $4$-GDD of type $2^{6}5^{1}11^{3}$ in Lemma 2.4
$\begin{array}[]{|l|l|l|l|l|}\hline\cr\\{c_{0},c_{2},c_{4},q_{0}\\}&\\{a_{0},a_{1},c_{2},s_{0}\\}&\\{a_{0},d_{4},f_{2},q_{0}\\}&\\{b_{0},c_{5},f_{1},p_{1}\\}&\\{c_{0},d_{2},f_{1},t_{2}\\}\\\
\hline\cr\\{e_{0},e_{2},e_{4},q_{0}\\}&\\{a_{0},a_{2},d_{1},t_{0}\\}&\\{a_{0},e_{5},f_{1},s_{3}\\}&\\{b_{0},c_{2},t_{1},t_{2}\\}&\\{c_{0},q_{1},s_{2},t_{1}\\}\\\
\hline\cr\\{s_{0},s_{2},s_{4},q_{0}\\}&\\{a_{0},b_{1},c_{5},s_{4}\\}&\\{a_{0},e_{2},q_{1},t_{2}\\}&\\{b_{0},d_{5},e_{4},t_{0}\\}&\\{d_{0},d_{1},p_{0},t_{3}\\}\\\
\hline\cr\\{f_{0},f_{1},t_{0},t_{3}\\}&\\{a_{0},b_{4},d_{2},s_{2}\\}&\\{a_{0},e_{4},t_{1},t_{5}\\}&\\{b_{0},e_{1},f_{2},s_{0}\\}&\\{d_{0},e_{1},e_{2},s_{4}\\}\\\
\hline\cr\\{f_{0},s_{0},s_{3},\infty\\}&\\{a_{0},b_{2},e_{1},p_{0}\\}&\\{b_{0},b_{1},d_{2},q_{1}\\}&\\{c_{0},c_{1},e_{5},t_{4}\\}&\\{d_{0},s_{2},s_{3},t_{4}\\}\\\
\hline\cr&\\{a_{0},b_{5},t_{3},\infty\\}&\\{b_{0},b_{2},s_{1},t_{5}\\}&\\{c_{0},d_{1},d_{5},s_{0}\\}&\\{e_{0},p_{0},s_{1},t_{4}\\}\\\
\hline\cr&\\{a_{0},c_{4},p_{1},s_{1}\\}&\\{b_{0},c_{1},e_{2},s_{2}\\}&\\{c_{0},d_{4},e_{2},\infty\\}&\\\
\hline\cr\end{array}$

######  Lemma 2.5.

There exists a $4$-GDD of type $2^{3}11^{4}$.

###### Proof.

The points in the $4$-GDD are

  * •

$a_{i},b_{i},c_{i},d_{i},e_{i},r_{i},t_{i}$ for $i\in{\mathbb{Z}}_{6}$;

  * •

$f_{i}$ for $i\in{\mathbb{Z}}_{3}$;

  * •

$p_{i},q_{i}$ for $i\in{\mathbb{Z}}_{2}$; and

  * •

$\infty$.

The groups are

  * •

$\\{a_{i},a_{i+3},b_{i},b_{i+3},c_{i},c_{i+3},d_{i},d_{i+3},e_{i},e_{i+3},f_{i}\\}$
and $\\{t_{i},t_{i+3}\\}$ for $i\in\\{0,1,2\\}$; and

  * •

$\\{p_{i},q_{i}\>:\>i\in{\mathbb{Z}}_{2}\\}\cup\\{r_{i}\>:\>i\in{\mathbb{Z}}_{6}\\}\cup\\{\infty\\}$.

Blocks are obtained by developing modulo $6$ the subscripts of the noninfinite
points in the base blocks given in Table 5. The first two blocks in the first
column generate two blocks each, while the third block in that column
generates one block. ∎

Table 5: Base blocks for the $4$-GDD of type $2^{3}11^{4}$ in Lemma 2.5
$\begin{array}[]{|l|l|l|l|l|l|}\hline\cr\\{a_{0},a_{2},a_{4},p_{0}\\}&\\{a_{0},a_{1},c_{2},r_{0}\\}&\\{a_{0},c_{4},f_{2},t_{2}\\}&\\{b_{0},b_{1},f_{2},r_{4}\\}&\\{c_{0},c_{1},r_{2},t_{0}\\}&\\{d_{0},f_{2},r_{3},t_{4}\\}\\\
\hline\cr\\{c_{0},c_{2},c_{4},q_{0}\\}&\\{a_{0},b_{1},d_{2},p_{1}\\}&\\{a_{0},d_{4},e_{2},r_{3}\\}&\\{b_{0},c_{4},e_{5},q_{1}\\}&\\{c_{0},d_{4},d_{5},t_{1}\\}&\\{e_{0},e_{1},p_{0},t_{1}\\}\\\
\hline\cr\\{f_{0},f_{1},f_{2},\infty\\}&\\{a_{0},b_{2},d_{1},q_{0}\\}&\\{a_{0},d_{5},q_{1},t_{5}\\}&\\{b_{0},c_{2},e_{4},r_{1}\\}&\\{c_{0},d_{2},e_{4},\infty\\}&\\{e_{0},f_{1},q_{1},t_{2}\\}\\\
\hline\cr&\\{a_{0},b_{5},e_{1},t_{0}\\}&\\{a_{0},e_{5},f_{1},r_{4}\\}&\\{b_{0},c_{1},p_{1},t_{3}\\}&\\{c_{0},d_{1},f_{2},p_{1}\\}&\\\
\hline\cr&\\{a_{0},b_{4},t_{4},\infty\\}&\\{a_{0},r_{1},t_{1},t_{3}\\}&\\{b_{0},c_{5},r_{5},t_{2}\\}&\\{d_{0},d_{2},r_{2},t_{1}\\}&\\\
\hline\cr&\\{a_{0},c_{5},e_{4},r_{2}\\}&\\{b_{0},b_{2},d_{4},r_{2}\\}&\\{b_{0},e_{1},t_{4},t_{5}\\}&\\{d_{0},e_{1},e_{5},r_{1}\\}&\\\
\hline\cr\end{array}$

##  3 Acknowledgements

This research used the computational cluster Katana supported by Research
Technology Services at UNSW Sydney. The third author acknowledges the support
from an Australian Government Research Training Program Scholarship and from
the School of Mathematics and Statistics, UNSW Sydney.

## ORCID

R. J. R. Abel: https://orcid.org/0000-0002-3632-9612
T. Britz: https://orcid.org/0000-0003-4891-3055
Y. A. Bunjamin: https://orcid.org/0000-0001-6849-2986
D. Combe: https://orcid.org/0000-0002-1055-3894

## References

  * [1] R.J.R. Abel, T. Britz, Y.A. Bunjamin and D. Combe, Group divisible designs with block size $4$ and group sizes $2$ and $5$, J. Combin. Des. 30 (2022), 367–383.
  * [2] R.J.R. Abel, T. Britz, Y.A. Bunjamin, D. Combe and T. Feng, Group divisible designs with block size 4 where the group sizes are congruent to $1\pmod{3}$, Discrete Math. 346 (2023), 113277.
  * [3] R.J.R. Abel, Y.A. Bunjamin and D. Combe, Some new group divisible designs with block size $4$ and two or three group sizes, J. Combin. Des. 28 (2020), 614–628.
  * [4] R.J.R. Abel, Y.A. Bunjamin and D. Combe, Existence of 4-GDDs with at most 50 points and $4$-GDDs of types $6^{s}3^{t}$ and $9^{s}3^{t}$, Discrete Math. 344 (2021), 112479.
  * [5] D. Kreher and D.R. Stinson, Small group divisible designs with block size four, J. Statist. Plann. Inference 58 (1997), 111–118.
  * [6] R. Rees and D.R. Stinson, Kirkman triple systems with maximum subsystems, Ars Combin. 25 (1988), 125–132.


[1]Sachini Karunarathne

These authors contributed equally to this work.

Sachini Karunarathne: methodology and writing, Kushani De Silva:
conceptualization, methodology, writing/editing, Sanjeewa Perera:
conceptualization

[1,2,3]Department of Mathematics, University of Colombo, Colombo, Sri Lanka

# Modeling and indexing the drought severity by augmenting SPI metric system
to incorporate multi-modal ground temperature data

<EMAIL_ADDRESS>Kushani De Silva<EMAIL_ADDRESS>Sanjeewa Perera<EMAIL_ADDRESS>*

###### Abstract

Drought is a global threat caused by the persistent challenges of climate
change. It is important to identify drought conditions based on the weather
variables and their patterns. In this study, we enhanced the Standardized
Precipitation Index (SPI) by integrating ground temperature data to develop a
more comprehensive metric for evaluating drought severity. Our metric offers a
dual assessment of drought severity, taking into account both the intensity of
the drought and its duration. We employ this evaluation in the primary paddy
cultivation region of Sri Lanka, with the aim of shedding light on the
prevailing drought conditions affecting paddy crops due to insufficient water
supply and prolonged periods of elevated temperatures. Additionally, we
calibrate our metric by aligning it with historical drought records and
subsequently compare the outcomes with those derived from the conventional SPI
index.

###### keywords:

Drought Severity, SPI, MSDI, Paddy, Copula, Multi-modal Distributions

##  1 Introduction

Droughts represent a substantial threat to nations and their inhabitants in
terms of economy, food security, health and much more. Currently, multiple
regions worldwide are suffering with droughts at varying degrees of severity.
The examination of droughts is of utmost importance, considering that the
likelihood of these events becoming more intense in the future is high,
primarily due to the persistent challenges posed by climate change. Concerning
food security and its relationship with climate change, rice, plays a pivotal
role. Rice serves as a staple food for more than half of the global population
[10], particularly in Asia, where it serves as a fundamental source of
calories and nutrition. Moreover, paddy fields contribute to environmental
sustainability by acting as carbon sinks, aiding in the mitigation of climate
change. In particular, paddy cultivation is of paramount importance to Sri
Lanka, both historically and in the present day. The country’s culture,
economy, social fabric and food security have deep ties to rice production [4,
29]. North Central province (NC) of Sri Lanka makes the largest contribution
to the country’s agriculture and enable 8.1% of country’s GDP and
approximately 107 kg per person of daily consumption necessity for Sri Lankan
population [13, 6, 8, 21, 18, 23]. Paddy cultivation is supported by the
weather of two main monsoons referred to as Yala and Maha - Maha seasons from
September to March and Yala season from May to August. NC experiences a
combination of both of these weather patterns at a given month due to its
geographical location.
Analyzing individual factors contributing to droughts falls short in offering
a comprehensive understanding of the complex ramifications of drought events
[35, 11]. Consequently, relying solely on metrics like Standardized
Precipitation Index (SPI) proves insufficient in gauging drought severity
because it primarily assesses the influence of rainfall [14, 28]. However,
temperature is a crucial parameter in discerning droughts since certain
temperature conditions can offset the effects of rainfall. Moreover, the
conventional SPI metric is based on the assumption that rainfall adheres to a
gamma distribution [19]. However, our research has unveiled that the gamma
distribution is insufficient for effectively modeling rainfall data. As a
result, we cast doubt on the appropriateness of employing this distribution
for the calculation of any drought index. This matter becomes particularly
pertinent in regions where meteorological variables do not conform to a single
distribution but exhibit multiple distributions due to fluctuating seasonal
weather patterns. We have encountered such a situation in the NC province,
where the influence of two distinct monsoon seasons has caused temperature
data to deviate from a unified probability distribution framework.
It is understood that hazards such as droughts are caused by extreme weather
conditions where conventional joint distributions would fail to model.
Therefore, it is important to choose a proper model that enables to model
extreme cases of weather that would pay more attention to the tails of the
respective marginal distributions. Copula is one of the robust frameworks that
enable modeling tail dependencies, and by that, it allows to capture the
consequences of extreme weathers that lie in tails of the distributions [15].
The selection of a suitable copula is based upon meticulous considerations
such as Kendall’s tau rank correlation measure and the shapes of the marginal
distributions of the selected variables [32]. The Kendall’s tau
$\left(\tau\right)$, measures the rank correlation providing insights into the
strength and direction of dependence between variables. Several examples of
modeling extreme cases of drought variables can be found at [24, 31] and
references therein. However, in these studies, variables such as drought
severity and drought duration have been used as the drought variables [32, 17,
16]. In contrast we build the index by directly using the weather data
avoiding the uncertainty in calculated drought variables.
In this study, we introduce an enhanced calculation for SPI that provides a
more precise representation of rainfall by utilizing a more fitting
probability distribution. Furthermore, we incorporate temperature data to
construct a more accurate metric for measuring drought. This involves
addressing the challenge of modeling ground temperature in regions
characterized by multiple weather seasons. Subsequently, we establish a
drought index capable of assessing drought severity from two perspectives, and
we compare our findings with those obtained using the traditional SPI index.
Our new metric’s validity is verified through an examination of historical
drought records, offering valuable insights into paddy-related drought
conditions.
The structure of this paper is as follows: Section 2 outlines the methodology
for fitting an appropriate copula to the marginal probability distributions of
selected weather variables, while Section 3 details the development of the
drought measuring metric and the process of its calibration followed by
Results and Discussion. The paper concludes with a summary in the final
section with valuable future avenues to explore.

##  2 Jointly dependent weather model

Driven by the rationale presented in the introduction, we have chosen NC as
our study area (latitude: $8.1996$ °N and longitude: $80.6327$ °E), covering a
span of 40 years from 1981 to 2021 to demonstrate the applicability of the
contributions of this work. The average rainfall and temperature conditions in
NC and favorable thresholds of these variables for the paddy are given in
Table 1. According to Table 1, although the favorable conditions of rainfall
and temperature agrees with the district wise monthly average values, NC has
experienced drought conditions numerous times in the past [13]. This
observation underscores the need for multivariate modeling, emphasizing the
importance of not relying solely on SPI and advocating for the inclusion of
prolonged periods in drought decision-making processes [26]. To achieve these
goals, we model the drought conditions using copula functions that enable
identifying joint structures between variables. The monthly average rainfall
(mm/day) and monthly temperature data (°C) were retrieved from the NASA
website [5]. In the next subsections, we demonstrate how the marginal
probability distributions are fitted to rainfall and temperature data
respectively to finally showcase the joint model using copula.

Table 1: Weather characteristics of the NC province’s districts are illustrated [7]. The favorable environmental conditions are also given for comparisons [9]. Although the average rainfall and temperature values lie in the favorable bounds, NC has faced droughts for paddy numerous times. The annual favorable values are not shown as they do not make any useful contribution to the study. District |  Monthly Average Rainfall (mm/day) |  Monthly Average Temperature (°C) |  Total Annual Rainfall |  Avg. Annual Temperature
---|---|---|---|---
Anuradhapura |  $3.6$ |  $28.3$ |  1533.4 mm |  28.2 ° C
Polonnaruwa |  $3.2$ |  $31.0$ |  1698.3 mm |  28.8 ° C
NC province |  $4.3$ |  $28.9$ |  $1589.0$ mm |  $27.3$ ° C
Favorable conditions |  $\left(2.74-4.11\right)$ |  $\left(21-27\right)$ |  |
Figure 1: Time series plots of monthly average temperature and rainfall in NC
for the period of 1981-2021. The favorable thresholds for paddy are marked
with dashed horizontal lines.

###  2.1 Estimating the marginal distributions

In order to determine the marginal probability distributions of each variable,
the rainfall and temperature data were fitted to different known probability
functions. Afterwards, the best fit was chosen based on AIC criteria
$\text{AIC}=2k-2\ln({\mathcal{L}}),$ where $k$ and ${\mathcal{L}}$ are
respectively the number of parameters and the likelihood function of the model
[2]. A likelihood function defined on the parameter space $\Omega$ is given
by,

$\mathcal{L}(\Theta\mid
x_{1},x_{2},\ldots,x_{n})=\prod_{i=1}^{n}f(x_{i};\Theta)$ (1)

where $\Theta$ is the vector of model parameters, $x_{i},i=1,\cdots,n$
represent the $n-$ sized dataset and $f(.)$ is the probability density
function of $x_{i}$. According to the AIC, the rainfall data (Fig. 2) agrees
the beta distribution rather than the gamma distribution (see Table 2). This
is evident that gamma distribution framework is lacking to represent rainfall.
Thus, this measurement has increased the accuracy of a conventional SPI value.
The mean of rainfall obtained from the fitted beta distribution is
$E(R)=0.128$ and it is outside of the favorable rainfall for paddy. Moreover,
the probability of NC rainfall (R) being favorable to paddy is $14.04$. i.e.
$p\left(R\in\left(2.74\text{ mm },4.11\text{ mm }\right)\right)=14.04\%$ which
indicates that 86% of the time, rainfall in NC is not favorable to paddy.

Table 2: Rainfall: fitted probability distributions and their corresponding AIC values for rainfall data. HG represents the Half-Gaussian distribution, i.e. $N(\mu,\sigma)$ where $\mu=0$ and $\sigma\in\left(0,\infty\right)$.  Probability distribution | Parameter estimates | AIC value
---|---|---
Beta $(\alpha_{R},\beta_{R})$ |  $\alpha_{R}=0.779$, $\beta_{R}=5.271$ | 2061.749
Gamma $(k_{R},\theta_{R})$ |  $k_{R}=0.861$, $\theta_{R}=3.690$ | 2104.167
HG $(\mu_{R},\sigma_{R})$ | $\mu_{R}=0,\sigma_{R}=4.821$ | 2265.897
Figure 2: Histogram of monthly average rainfall and temperature data. Rainfall
data are rightly skewed distributed in the range of $\left(0-18.41\right)$
mm/day. The yellow line shows the best-fitted beta distribution for rainfall.
Histogram of monthly average temperature shows a multi-modal behavior where a
local peak within the favorable temperatures (vertical dashed lines) and the
global peak is outside the favorable temperature. The data are distributed in
the range of (22.98-33.74) °C. The yellow solid line shows the GMM fitted to
the temperature data.

The histogram of temperature data depicted in Fig. 2 shows clear evidence of
non-Gaussian multi-modals. The two peaks of the histogram fairly agrees to
that of a Gaussian, i.e. assuming data are generated by a combination of two
Gaussian distributions. The Gaussian Mixed Model (GMM) of two components is
fitted to the temperature data and it is mathematically denoted as below:

$f(x|\Theta_{i})=\sum_{i=1}^{n}w_{i}\cdot\mathcal{N}(\mu_{i},\sigma_{i}^{2}),$
(2)

where $\Theta_{i}=\left(w_{i},\mu_{i},\sigma_{i}\right)$. Here,
$w_{i},\sum_{i}w_{i}=1$ denotes the weight assigned to each component
distributed with mean $\mu_{i}$ and $\sigma_{i}$. The parameters were
estimated using the iterative Expectation-Maximization (EM) algorithm where
initial guesses for the parameters were obtained by K-clustering algorithm.
The parameter estimates for $\Theta_{i}=\left(w_{i},\mu_{i},\sigma_{i}\right)$
were estimated for two components respectively as
$\Theta_{1}=\left(0.31,24.8586,0.8015\right)$ and
$\Theta_{2}=\left(0.69,28.6589,1.3191\right)$. The final fitted density is
given in Fig. 2.

###  2.2 Copula fitting

A copula is a function that joints several one-dimensional distribution
functions to a multivariate distribution function by capturing their
dependency structure [33]. For a bivarite case with variables $X$, $Y$ and
marginal distribution functions $F_{X}(x)$, $F_{Y}(y)$ then there exists a
copula $C$ such that,

$P(X\leq x,Y\leq y)=C(F_{X}(x),F_{Y}(y)).$ (3)

The copula $C$ is fitted to the transformed uniform data ($U_{X}$ and $U_{Y}$)
of the marginals. PIT (see Theorem 1) gurantees the uniformity of the random
variables found by transforming marginal distribution functions $F_{X}(x)$,
$F_{Y}(y)$ [22].

######  Theorem 1.

Probability Integral Transform Theorem (PIT)
If $X$ is a continuous random variable with distribution function $F$, then
the random variable $U=F(X)$ (the probability integral transform of $X$) is
uniformly distributed on the unit interval $I=[0,1]$; or equivalently,
$P[F(X)\leq t]=t$ for all $t\in I$.

The marginal densities obtained in the previous section is then transformed
into uniform distributions using PIT to identify a suitable copula function.
Obtaining a best-fit copula is done in two steps: (a) finding the appropriate
copula from set of pre-defined copula functions using a goodness of fit test,
(b) optimizing the tune-in parameter of the selected copula. The value of
Kendall’s $\tau\in\left[-1,1\right]$ was $-0.00711$ for the uniform data of
rainfall ($u$) and temperature ($v$). Based on this value of $\tau$, Frank and
Farlie-Gumbel Morgenstern (FGM) copulas were suitable as they are capable of
capturing low dependencies. The selected copulas, their distributions and
density functions are given in Table 3 [25, 32]. To find the proper dependency
structure between weather axes, these two copulas were fitted to uniform data.

Table 3: Set of appropriate copula functions according to the resulting low dependency Kendall’s tau rank correlation. The distribution function and density function of each copula is given by uppercase and lowercase $c(.)$ respectively. The copula parameter is represented by $\theta$ and uniform data for rainfall and temperature are given by $u,v$ respectively. The estimated copula parameter and the p-value of the goodness of fit values ($S_{n}$) are also given. | Frank | Farlie-Gumbel-Morgenstern (FGM)
---|---|---
$C(u,v)$ | $\dfrac{1}{\theta}\left(1+\dfrac{(e^{-\theta u}-1)(e^{-\theta v}-1)}{e^{-\theta}-1}\right)$ | $uv\left(1+\theta(1-u)(1-v)\right)$
$c(u,v)$ | $-\dfrac{\theta e^{-\theta(u+v)}(e^{-\theta}-1)}{(e^{-\theta(u+v)}-e^{-\theta u}-e^{-\theta v}+e^{-\theta})^{2}}$ | $1+\theta(1-2u)(1-2v)$
Estimated $\theta$ | $-0.0013$ | $-0.03195$
p-value of $S_{n}$ | $0.6389$ | $0.6968$

The most appropriate copula is selected by following the hypothesis testing
given below:

$H_{0}$

: The given copula $C_{0}$ belongs to the chosen parametric copula family $C$.

$H_{1}$

: The given copula $C_{0}$ does not belong to the chosen parametric copula
family $C$.

The test-statistic to select the correct hypothesis is the Cramer-von Mises
statistic $S_{n}$ for $n$ number of data points [12, 27],

$S_{n}=\sum_{i=1}^{n}(C_{n}({u}_{in},{v}_{in})-C_{\theta_{n}}({u}_{in},{v}_{in}))^{2},\quad$
(4)

where $C_{n}=\frac{1}{n}\sum_{i=1}^{n}1(u_{in}\leq u,v_{in}\leq v)$, and
$u_{in},v_{in}$ are given by $\frac{1}{n+1}\sum_{j=1}^{n}1(X_{j}\leq X_{i})$
and $\frac{1}{n+1}\sum_{j=1}^{n}1(Y_{j}\leq Y_{i})$ respectively. The
empirical copula $C_{n}$ was obtained from the pseudo-observations
$\left(u_{in},v_{in}\right)_{i=1}^{n}$ calculated from
$(X_{1},Y_{1}),...(X_{n},Y_{n})$ where $X$ and $Y$ are the rainfall and
temperature data. According to Table 3, the p-values of both copula are
significant at 5% level. The strongest copula is therefore selected to be FGM
as it has the highest p-value [27]. The FGM copula density plot and contour
plot is given in Fig. 3.

Figure 3:  (a) Fitted FGM copula density and (b) Fitted FGM copula
distribution at $\theta=-0.0319$ . The axis labels $u,v$ respectively show the
transformed rainfall and temperature data.

##  3 Developing the drought severity metric

In this section we showcase how the copula model is used to calculate the
severity metric. The joint distribution with cumulative joint probability $p$
for weather data can be expressed with the FGM copula (see Eq. (5)).

$P(X\leq x,Y\leq y)=C[F(X),G(Y)]=p,$ (5)

where $F(X)$ and $G(Y)$ are the marginal cumulative distribution functions of
rainfall ($X$) and temperature ($Y$). In [14], they have introduced a method
of calculating their MSDI metric using,

$MSDI=\psi^{-1}(p),$ (6)

where $\psi$ is the standard normal distribution function. We simply borrow
Eq. (6) to calculate MSDI for our own variables. The metric calculated from
Eq. (6) gurantees to generate values in the range of
$\left(-\infty,\infty\right)$ by the defined range of standard normal
distribution. The infinity for the index is approximately taken as 5 due to
the fact that 99% of the probability of the normal distribution is between
$\left(-4,4\right)$. Negative and positive MSDI values respectively indicate
dry climate conditions (drought) and wet climate conditions. MSDI near zero
signify the normal climate conditions [14]. The indices generated in the
previous section were then re-computed for different time periods to showcase
the prolonged period of certain weather dependencies. Thus the two dimensions
of the drought assessment system are MSDI index and the prolonged period. Four
significant time periods were identified for T as 3-month, 6-month, 9-month
and 12-month. When the prolonged period, T is higher, the dry conditions from
MSDI reduced because less dry conditions for a longer period can cause for
severe droughts. The n-month MSDI is computed as,

$n\text{-month MSDI}=\frac{1}{n}\sum_{k=1}^{n}\text{1-month MSDI}_{k}.$ (7)

###  3.1 System calibration

The MSDI values generated for different T periods are categorized based on
identified thresholds to identify normal conditions and two different dry
conditions - namely dry and extremely dry. The thresholds identified can be
found at Table 4. For the two dry conditions, the thresholds of MSDI in
different T periods are depicted in Fig. 4 for more clarity. According to the
Fig. 4, when the prolonged period, T increases, the severity of dry weather
conditions decreases indicating severity $\propto\left(T\times MSDI\right)$.
Moreover, when comparing the dry conditions vs extremely dry conditions, the
levels of thresholds in MSDI are smaller for extremely dry conditions. The
threshold to differentiate dry from extremely dry increases as T increases. In
summary these trends indicate that we can reliably identify a drought
situation when there is a combination of high temperatures and low rainfall
persisting over an extended period. However, it’s important to recognize that
when we examine shorter time spans, it becomes more challenging to swiftly
determine whether a drought condition is present or not. Therefore, the
duration of the observed climate conditions significantly influences our
ability to classify them as a drought.

Figure 4: The trend of thresholds of normal, dry and extremely dry conditions along with prolonged period is depicted.  Table 4: Two-dimensional threshold values for MSDI to identify drought conditions are given.  Drought conditions | MSDI range
---|---
3-month | 6-month | 9-month | 12-month
Extremely dry | $\left(-\infty,-3.2\right)$ | $\left(-\infty,-2.5\right)$ | $\left(-\infty,-2.1\right)$ | $\left(-\infty,-1.7\right)$
Dry | $\left(-3.2,-1.3\right)$ | $\left(-2.5,-0.8\right)$ | $\left(-2.1,-0.4\right)$ | $\left(-1.7,-0.3\right)$
Normal | $\left(-1.3,+\infty\right)$ | $\left(-0.8,+\infty\right)$ | $\left(-0.4,+\infty\right)$ | $\left(-0.3,+\infty\right)$

##  4 Results and discussion

In this section, we discuss the applicability of the new metric developed in
identifying the severity of the paddy-related drought conditions and to
compare them against SPI thresholds. The new metric, MSDI is depicted in Fig.
5 whereas the SPI is depicted in Fig. 6. Fitting this long-term precipitation
record to a gamma probability distribution is the first step. Once the gamma
probability distribution is determined, the cumulative probability of observed
precipitation is computed and then inverse transformed by a standard normal
distribution with mean $0$ and variance $1$. The resulting index is the SPI
[20]. In Fig. 5, if an MSDI value extend beyond red (bold) line, it is
indicative of an extreme dry event whereas if the index value falls between
the red (thin) line and red (bold) line it is categorized as a dry condition.
The rest can be categorized under the normal/healthy conditions (see Table 4).

Sri Lanka experienced a widespread drought in 1992 with over 18 drought
incidents recorded [34]. According to Fig. 5, this drought is shown as
extremely dry in 3, 6, 9 and 12 months time scales but not very significant in
9 or 12 months. This confirms that the drought of 1992 was characterized by
numerous minor drought events recorded in sub-intervals, rather than a single
continuous drought episode. However, the SPI failed to identify this drought
episode in its 6, 9, or 12 months SPI charts in Fig. 6. The CBSL annual report
in [3] states that the paddy production is declined in the cultivation year
1982 due to adverse weather conditions. Additionally, it states that it’s
worth noting that there is a significant 13% drop in paddy production during
the Yala season of 1983 when compared to the previous Yala season. As per this
report, this notable decline in paddy production during Yala $1983$ shows the
abandoned or suffering crop damage due to the drought conditions that
prevailed during the first half of the year [3]. Furthermore, based on another
historical record, the period from 1982-1983 is categorized as a dry period
[34]. Our developed MSDI, computed for 3, 6, 9, and 12 months, consistently
reflects these dry conditions for $1982-1983$ period, as it consistently
hovers above the threshold of extreme dryness (represented by the red (bold)
line) and below the dry limit (red (thin) line). The length of the peaks of
$1983$ is longer than $1982$ as in Fig. 5, and it aligns with the statement of
the CBSL annual report-1983. That is, compared to $1982$ Yala season, there’s
a higher paddy loss in $1983$. Comparing the SPI, in 3 months it identifies
$1983$ as an extreme drought situation while $1982$ as a dry situation.
Furthermore, in $6,9$ timeframes SPI only denoted $1983$ as a dry situation.
In long run, SPI fails in identifying the dry scenarios for both timeframes as
indicated in Fig. 6.

Referring the historical recordings of CBSL annual report -1987, it emphasizes
that the decline of the paddy production in $1987$ was mainly due to the
prolonged drought conditions during the cultivation year 1986-1987 [3].
Further, it highlights that this was the second consecutive year in which
paddy production dropped where the estimated 1.39 million metric tons was a
substantially declined by 18%. If we consider the developed MSDI (given in
Fig. 5), in long run basically for 12 months period it clearly depicts a dry
scenario for the period of 1986-1987 where as SPI shows a completely normal
situation in long run as its peaks are above its standard thresholds. In both
6, 9 months cases MSDI was successful in identifying 1987 as a drought
condition which proves the long term existence of the drought in this
cultivation period. SPI is only succeeded in identifying the drought scenario
in the short term that is in 3 months case. Over the long run, the SPI
consistently falls outside its standard threshold limits, indicating its
ineffectiveness in identifying drought situations in prolong period as
indicated in Fig. 6.

Due to the drought condition prevailed in the 1988-1989 cultivation era, the
paddy production has been significantly reduced according to the CBSL-1989
report [3]. Moreover, it highlights that, fall in paddy production was
reflected in both 1988-1989 Maha and 1989 Yala, seasons. Considering our MSDI
for 9 months period it identifies a dry period for the cultivation year
1988-1989, but SPI is not capable in identifying the drought scenario in this
timeframe as its peaks lies outside its standard threshold limit, but SPI
identify this drought event in 3 months as shown in Fig. 6.

For the period of 1995-1996 CBSL annual report -1996 in [3] states that the
paddy production in Anuradhapura district dropped by 61% from 190,000 metric
tons produced during Maha 1994-1995 to 74,000 metric tons in Maha 1995-1996
due to the occurred drought situation. Referring the MSDI for the cultivation
year 1996 it is clearly visible that all the peaks lie below the extremely dry
situation for 3, 6, 9, 12 timeframes. This depicts the existence of this
drought event for a prolonged period.This suggests that the drought’s severity
is significantly greater. This observation aligns with the CBSL report, which
notes that paddy production, reaching its peak in 1995, experienced a sharp
decline of 27% to reach 2.1 million metric tons—the lowest recorded level
since 1979 [3]. However, SPI categorizes this drought occurrence as a dry
scenario. Consequently, when comparing the two indices, MSDI emerges as the
more dependable tool for recognizing this drought situation in accordance with
historical records.

According to the CBSL-2001 report, the extents sown and harvested paddy during
both Maha, Yala seasons dropped due to adverse weather conditions that
prevailed during the year [3]. The MSDI calculations for the 3, 6, and 9
months periods effectively pinpoint a dry condition, and over an extended
duration, the MSDI surpasses the extremely dry threshold, indicating the
persistence of this drought condition for a prolonged period (see Fig. 5).
However, it’s worth noting that SPI does not accurately represent this dry
condition for this cultivation year as shown in Fig. 6. Similarly in the
cultivation year 2004, the MSDI exceeds the extremely dry threshold for all
timeframes. However this can be proven using the historical record, CBSL-2004
annual report [3]. The report indicates that the agriculture sector’s
contribution to the GDP declined from 19% in 2003 to 18% in 2004, primarily
due to reduced paddy production. This decline is attributed to drought
conditions in the Anuradhapura and Kurunegala districts during the Maha
season, resulting from the failure of the Northeast monsoon and the delayed
onset of the Southwest monsoon in the Yala season [3].

For year 2014, the CBSL 2014 annual report in [3] highlights a slowdown in the
agriculture sector, with a growth rate of only 0.3% compared to 4.7% in the
previous year. This deceleration is primarily attributed to adverse weather
conditions that prevailed throughout the year. Notably, while the SPI have not
indicated this adverse drought scenario accurately, the MSDI effectively
denotes peaks for drought conditions for 3, 6, 9, and 12 months periods.
Additionally, as documented in historical records, this drought condition
persisted for an entire year, as evidenced by the visible peaks even in the 12
months MSDI panel (Fig. 5). During the drought period spanning 2016-2017, the
MSDI peaks fell within the range between the extremely dry and dry threshold
limits. The CBSL 2016 annual report in [3] highlights the significant impact
of adverse weather conditions, particularly the prolonged drought, on
agricultural performance throughout that year. Interestingly, while the MSDI
accurately identifies this drought situation (see Fig. 5), the SPI, in
contrast, does not recognize it and instead categorizes it as a normal
scenario (see Fig. 6). All in all, the calculated two-dimensional MSDI
incorporating rainfall with temperature outperforms the SPI when calibrated
against historical drought records.

Figure 5: The MSDI for 3, 6, 9, 12 months for the period of 1981-2021 are
shown respectively in (a), (b), (c) and (d). The highest downward peaks
represent the years of the occurrence of severe droughts. The red(bold) and
red(thin) dash lines represent threshold limits for extremely dry and dry
conditions respectively. Figure 6: The MSDI and SPI for 3, 6, 9, 12 months for
the period of 1981-2021 are shown respectively in (a), (b), s(c) and (d). The
highest downward peaks represent years of the occurrence of severe droughts.
The red and green dash lines represent threshold limits of extremely dry
conditions of MSDI and SPI respectively.

##  5 Conclusions

As outlined in the introduction, droughts result from two primary factors: (1)
prolonged periods of elevated temperatures and (2) extreme precipitation
events. Consequently, our drought severity assessment system takes into
account the effects of these weather variables over extended time periods. The
development of this two-dimensional index system involved calibration using
historical drought events, demonstrating that the index aligns with the
severity experienced during those drought occurrences.
The North Central (NC) region’s semi-arid climate makes it highly susceptible
to droughts, which pose significant challenges to the agricultural sector.
Droughts occur during prolonged periods of reduced precipitation, aggravated
by elevated temperatures. As a result, the paddy cultivation industry in this
area is heavily affected. NC has a history of experiencing numerous droughts,
including instances in 1935-1937, 1947-1949, 1953-1956, 1965, 1974-1977,
1982-1983, 1986-1987, 1992, 1995-1996, 2001, 2004, 2013-2014, and 2016-2017
[13, 30, 1]. In this study, we employ a weather data modeling approach to
assess drought severity, utilizing the FGM copula due to its capacity to
capture extreme weather events through tail dependencies. Our analysis focuses
on two key factors contributing to droughts, namely rainfall and temperature,
over a 40-year period spanning from 1981 to 2021. The marginal distributions
derived for rainfall and temperature are represented respectively by the beta
distribution (Fig. 2) and a Gaussian mixed model (Fig. 2). Following the
fitting of the copula to the transformed data, we compute the MSDI. These MSDI
values serve as indicators of drought severity, and we establish thresholds
for different levels of seventies, including normal, dry, and extremely dry
conditions (see Table 4).
Furthermore, the developed two-dimensional index system has been compared with
the already existing, most widely used drought monitoring index SPI [14]. It’s
noteworthy that the SPI is built upon the analysis of precipitation data only
[28]. The common approach of developing the SPI involves fitting a gamma
distribution to the precipitation data [19]. However, in our study we
integrate both precipitation and ground temperature to calculate MSDI. While
constructing the SPI, we followed the typical procedure, that is, fitting a
gamma distribution to the rainfall data though it best fits for a beta
distribution as detailed earlier. Though, for MSDI we employed a beta
distribution for rainfall according to the fitting choice of AIC value.
Afterwards, evaluating the developed MSDI against SPI, the historical data
played a crucial role in comparing and decision-making. Upon comparing these
two indices MSDI consistently outperforms SPI. Furthermore, we can conclude
that MSDI works more accurately for longer time frames.
In conclusion, we emphasize the significance of incorporating temperature data
into the index calculation, as opposed to relying solely on rainfall, as done
in the SPI. We illustrate the superior performance of MSDI compared to SPI.
Furthermore, we underscore the importance of accurately characterizing the
marginal distributions of weather variables. Our MSDI system’s validity is
affirmed through its successful validation against historical drought events.
When contrasting with existing drought severity evaluations found in the
literature, these assessments typically employ drought variables like drought
duration and drought intensity to construct frequency distributions. For
prospective research endeavors, the novel MSDI can be employed to compute
these drought variables, forming frequency distributions for drought
occurrences. This approach allows for a comparison of accuracy against
conventional models.

## References

  *   * Abeysingha and Rajapaksha [2020] Abeysingha N, Rajapaksha U (2020) Spi-based spatiotemporal drought over sri lanka. Advances in Meteorology 2020:10
  * Akaike [1974] Akaike H (1974) A new look at the statistical model identification. IEEE Transactions on Automatic Control 19(6):716–723
  * CBSL [2021a] CBSL (2021a) Cbsl annual report. Retrieved at: https://www.cbsl.gov.lk/en/publications/economic-and-financial-reports/annual-reports
  * CBSL [2021b] CBSL (2021b) Central bank of sri lanka annual report 2021. Retrieved at: https://www.cbsl.gov.lk/en/publications/economic-and-financial-reports/annual-reports/annual-report-2021
  * Center [2021] Center NLR (2021) Nasa power data access viewer. Retrieved at: https://power.larc.nasa.gov/data-access-viewer/
  * Chandrasiri et al [2020] Chandrasiri S, Galagedara L, Mowjood M (2020) Impacts of rainfall variability on paddy production: A case from bayawa minor irrigation tank in sri lanka. Paddy and Water Environment 18:443–454
  * DCS [2023] DCS (2023) Statistical pocket handbook: Department of census and statistics. Retrieved at: http://www.statistics.gov.lk/ref/PocketBook2022_Si_Ta
  * DMC [2017] DMC (2017) Disaster information management system-sri lanka. Retrieved at: http://www.desinventar.lk/
  * EOS [2023] EOS (2023) Growing rice: Sowing, cultivating, and harvesting. Retrieved at: https://eos.com/blog/how-to-grow-rice
  * Fukagawa and Ziska [2019] Fukagawa NK, Ziska LH (2019) Rice: Importance for global nutrition. Journal of nutritional science and vitaminology 65(Supplement):S2–S3
  * Ganguli and Reddy [2012] Ganguli P, Reddy M (2012) Risk assessment of droughts in gujarat using bivariate copulas. Water resources management 26:3301–3327
  * Genest et al [2009] Genest C, Rémillard B, Beaudoin D (2009) Goodness-of-fit tests for copulas: A review and a power study. Insurance: Mathematics and Economics 44:119–213
  * Gunawardhana and Dharmasiri [2016] Gunawardhana L, Dharmasiri L (2016) Drought hazard and managing its impacts through the disaster management approach: A study in the north central province of sri lanka. Paper presented at the International Research Symposium Rajarata University of Sri Lanka
  * Hao and AghaKouchak [2014] Hao Z, AghaKouchak A (2014) Multivariate standardized drought index: a parametric multi-index model. Advances in Water Resources 57:12–18
  * Haugh [2016] Haugh M (2016) An introduction to copulas. IEOR E4602: quantitative risk management Lecture notes Columbia University
  * Khan et al [2021] Khan MA, Faisal M, Hashmi M, et al (2021) Modeling drought duration and severity using two-dimensional copula. Journal of Atmospheric and Solar-Terrestrial Physics 214:105530
  * Lee et al [2013] Lee T, Modarres R, Ouarda T (2013) Data-based analysis of bivariate copula tail dependence for drought duration and severity. Hydrological Processes 27(10):1454–1463
  * Manesha et al [2015] Manesha S, Vimukthini S, Premalal K (2015) Develop drought monitoring in sri lanka using standard precipitation index (spi). Sri Lanka J Meteorol 1:64–71
  * McKee et al [1993a] McKee T, Doesken N, Kleist J (1993a) The relationship of drought frequency and duration to time scales. Paper presented at the Eighth Conference on Applied Climatology, 17-22 January 1993, Anaheim, California
  * McKee et al [1993b] McKee TB, Doesken NJ, Kleist J, et al (1993b) The relationship of drought frequency and duration to time scales. In: Proceedings of the 8th Conference on Applied Climatology, California, pp 179–183
  * MoE [2021] MoE (2021) Ministry of environment-sri lanka. Retrieved at: https://env.gov.lk/web/index.php/en/
  * Nelsen [2006] Nelsen RB (ed) (2006) An introduction to copulas. Springer Science & Business Media, New York
  * Nianthi [2019] Nianthi R (2019) Farmers’ responses to drought: Dry zone of sri lanka:(case study in medirigiriya). Case Studies Journal 5
  * Otkur et al [2021] Otkur A, Wu D, Zheng Y, et al (2021) Copula-based drought monitoring and assessment according to zonal and meridional temperature gradients. Atmosphere 12(8):1066
  * Pandey et al [2018] Pandey P, Das L, Jhajharia D, et al (2018) Modelling of interdependence between rainfall and temperature using copula. Modeling Earth Systems and Environment 4:867–879
  * Papadopoulos et al [2021] Papadopoulos C, Spiliotis M, Gkiougkis I, et al (2021) Fuzzy linear regression analysis for groundwater response to meteorological drought in the aquifer system of xanthi plain, ne greece. Journal of Hydroinformatics 23(5):1112–1129
  * Reddy and Ganguli [2012] Reddy M, Ganguli P (2012) Bivariate flood frequency analysis of upper godavari river flows using archimedean copulas. Water Resources Management 26(14):3995–4018
  * Reyes et al [2022] Reyes L, Rangel H, Herazo L (2022) Adjustment of the standardized precipitation index (spi) for the evaluation of drought in the arroyo pechelín basin, colombia, under zero monthly precipitation conditions. Atmosphere 13(2):236
  * Senanayake and Premaratne [2016] Senanayake S, Premaratne S (2016) An analysis of the paddy/rice value chains in sri lanka. Asia-Pacific Journal of Rural Development 26(1):105–126
  * Senatilleke et al [2023] Senatilleke U, Sirisena J, Gunathilake M, et al (2023) Monitoring the meteorological and hydrological droughts in the largest river basin (mahaweli river) in sri lanka. Climate 11(3):57
  * Seyedabadi et al [2020] Seyedabadi M, Kavianpour M, Moazami S (2020) Multivariate drought risk analysis based on copula functions: a case study. Water Supply 20(6):2375–2388
  * Shiau [2006] Shiau J (2006) Fitting drought duration and severity with two-dimensional copulas. Water resources management 20(795–815):795–815
  * Sklar [1959] Sklar A (1959) Fonctions de r’epartition ‘a n dimensions et leurs marges. Publications de l’Institut de Statistique de l’Universit’e de Paris 8:229–231
  * Tuong and Bouman [2003] Tuong T, Bouman B (2003) Rice production in water-scarce environments. Water productivity in agriculture: Limits and opportunities for improvement 1:13–42
  * Wang et al [2019] Wang F, Wang Z, Yang H, et al (2019) Copula-based drought analysis using standardized precipitation evapotranspiration index: A case study in the yellow river basin, china. Water 11(6):1298