Nuclear Safety
NEA/CSNI/R(2021)10
June 2024
www.oecd-nea.org

Benchmark on External Events

Hazard Frequency and Magnitude
Statistical Modelling

Nuclear Energy Agency

NEA/CSNI/R(2021)10
Unclassified

English text only
11 June 2024

NUCLEAR ENERGY AGENCYCOMMITTEE ON THE SAFETY OF NUCLEAR INSTALLATIONSBenchmark on External Events Hazard Frequency and Magnitude Statistical
Modelling

This document is available in PDF format only.

JT03545830
OFDE

This document, as well as any data and map included herein, are without prejudice to the status of or sovereignty over any territory,
to the delimitation of international frontiers and boundaries and to the name of any territory, city or area.

2 | NEA/CSNI/R(2021)10
ORGANISATION FOR ECONOMIC CO-OPERATION AND DEVELOPMENTThe OECD is a unique forum where the governments of 38 democracies work together to address the economic,
social and environmental challenges of globalisation. The OECD is also at the forefront of efforts to understandand to help governments respond to new developments and concerns, such as corporate governance, theinformation economy and the challenges of an ageing population. The Organisation provides a setting where
governments can compare policy experiences, seek answers to common problems, identify good practice and work
to co-ordinate domestic and international policies.
The OECD member countries are: Australia, Austria, Belgium, Canada, Chile, Colombia, Costa Rica,
Czechia, Denmark, Estonia, Finland, France, Germany, Greece, Hungary, Iceland, Ireland, Israel, Italy, Japan,
Korea, Latvia, Lithuania, Luxembourg, Mexico, the Netherlands, New Zealand, Norway, Poland, Portugal, theSlovak Republic, Slovenia, Spain, Sweden, Switzerland, Türkiye, the United Kingdom and the United States. The
European Commission takes part in the work of the OECD.
OECD Publishing disseminates widely the results of the Organisation’s statistics gathering and research on
economic, social and environmental issues, as well as the conventions, guidelines and standards agreed by its
members.

NUCLEAR ENERGY AGENCYThe OECD Nuclear Energy Agency (NEA) was established on 1 February 1958. Current NEA membershipconsists of 34 countries: Argentina, Australia, Austria, Belgium, Bulgaria, Canada, Czechia, Denmark, Finland,
France, Germany, Greece, Hungary, Iceland, Ireland, Italy, Japan, Korea, Luxembourg, Mexico, the Netherlands,
Norway, Poland, Portugal, Romania, Russia (suspended), the Slovak Republic, Slovenia, Spain, Sweden,
Switzerland, Türkiye, the United Kingdom and the United States. The European Commission and the International
Atomic Energy Agency also take part in the work of the Agency.
The mission of the NEA is:– to assist its member countries in maintaining and further developing, through international co-operation, thescientific, technological and legal bases required for a safe, environmentally sound and economical use of
nuclear energy for peaceful purposes;– to provide authoritative assessments and to forge common understandings on key issues as input togovernment decisions on nuclear energy policy and to broader OECD analyses in areas such as energy andthe sustainable development of low-carbon economies.
Specific areas of competence of the NEA include the safety and regulation of nuclear activities, radioactive
waste management and decommissioning, radiological protection, nuclear science, economic and technicalanalyses of the nuclear fuel cycle, nuclear law and liability, and public information. The NEA Data Bank provides
nuclear data and computer program services for participating countries.
This document, as well as any data and map included herein, are without prejudice to the status of or sovereignty over any territory, to thedelimitation of international frontiers and boundaries and to the name of any territory, city or area.
Corrigenda to OECD publications may be found online at: www.oecd.org/about/publishing/corrigenda.htm.
© OECD 2024

You can copy, download or print OECD content for your own use, and you can include excerpts from OECD publications, databases and multimedia products in yourown documents, presentations, blogs, websites and teaching materials, provided that suitable acknowledgement of the OECD as source and copyright owner is given.
All requests for public or commercial use and translation rights should be submitted to neapub@oecd-nea.org. Requests for permission to photocopy portions of thismaterial for public or commercial use shall be addressed directly to the Copyright Clearance Center (CCC) at info@copyright.com or the Centre français d'exploitation
du droit de copie (CFC) contact@cfcopies.com.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 3

Committee on the Safety of Nuclear Installations (CSNI)The Committee on the Safety of Nuclear Installations (CSNI) addresses NEA
programmes and activities that support maintaining and advancing the scientific andtechnical knowledge base of the safety of nuclear installations.
The Committee constitutes a forum for the exchange of technical information and forcollaboration between organisations, which can contribute, from their respective
backgrounds in research, development and engineering, to its activities. It has regard to
the exchange of information between member countries and safety R&D programmes ofvarious sizes in order to keep all member countries involved in and abreast ofdevelopments in technical safety matters.
The Committee reviews the state of knowledge on important topics of nuclear safety
science and techniques and of safety assessments, and ensures that operating experienceis appropriately accounted for in its activities. It initiates and conducts programmes
identified by these reviews and assessments in order to confirm safety, overcome
discrepancies, develop improvements and reach consensus on technical issues ofcommon interest. It promotes the co-ordination of work in different member countries
that serve to maintain and enhance competence in nuclear safety matters, including theestablishment of joint undertakings (e.g. joint research and data projects), and assists inthe feedback of the results to participating organisations. The Committee ensures thatvaluable end-products of the technical reviews and analyses are provided to members ina timely manner, and made publicly available when appropriate, to support broader
nuclear safety.
The Committee focuses primarily on the safety aspects of existing power reactors, othernuclear installations and new power reactors; it also considers the safety implications ofscientific and technical developments of future reactor technologies and designs.
Further, the scope for the Committee includes human and organisational research
activities and technical developments that affect nuclear safety.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

4 | NEA/CSNI/R(2021)10

Foreword
The modelling of external hazards encompasses different technical aspects, dependingon the type of hazard, though a general feature of hazards is that they produce off-normal
conditions that can impact nuclear installations. There is also a coupling of the hazard
with its associated risk analysis models. A risk analysis contains a set of scenarios,
frequencies and associated consequences, developed in such a way as to inform
decisions. A scenario contains an initiating event and (usually) one or more subsequent
events leading to an end state that reflects the issue of concern. The objective of thisbenchmark study is to focus on the initiating event by facilitating an exercise on thestatistical modelling for assessing hazard frequency and magnitude for external events
risk assessment. This benchmark study report provides details (data and overall
objectives) for the benchmarking exercise by specifying synthetic data for a hypotheticalexternal event (e.g. precipitation, extreme temperatures and high winds). The analysis
steps and modelling results of the benchmark participants are provided and summarised.
Overall conclusions from these submissions are described to gain insights from theactivity.
This report was approved by the Nuclear Energy Agency (NEA) Committee on the
Safety of Nuclear Installations (CSNI) at the 69th meeting of the CSNI held on 2-3 June
2021 (NEA/SEN/SIN(2021)1, not publicly available).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 5

Acknowledgements
Special thanks are extended to the following individuals for their contributions to thisbenchmarking exercise – Olle Räty and Marko Laine (Finnish Meteorological Institute);
Herve Cordier, Anne-Hélène Dutfoy Lebrun, Jeanne Ady Demgne and Luc Guinard(Électricité de France); Vincent Rebour, Yasser Hamdi and Yann Richet (Institut deRadioprotection et de Sûreté Nucléaire, France); Minkyu Kim, Beomjin Kim and KunYeon Han (Korea Atomic Energy Research Institute); Stef Carelsen (Autoriteit
Nucleaire Veiligheid en Stralingsbescherming, Netherlands); and Curtis Smith, Kurt
Vedros, Zhegang Ma, Jing Zhang and Andrea Mack (Idaho National Laboratory, UnitedStates).
Thanks are also due to all the participants who submitted their analysis for considerationand the writers and reviewers of this document.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

6 | NEA/CSNI/R(2021)10

Table of contents
List of abbreviations and acronyms.......................................................................................... 12
Executive summary .................................................................................................................... 13
1.

Introduction ......................................................................................................................... 15
1.1
1.2
1.3
1.4

2.

Benchmark problem description ....................................................................................... 19
2.1
2.2
2.3

3.

Background .............................................................................................................. 15
Objective .................................................................................................................. 16
Overview of the exercises ........................................................................................ 17
Glossary .................................................................................................................... 18
Data generation via synthetic models ....................................................................... 19
Case 1 – Known model producing the synthetic data .............................................. 19
Case 2 – Unknown model producing the synthetic data .......................................... 21

Benchmark outcomes .......................................................................................................... 23
3.1
Benchmark completion............................................................................................. 23
3.2
Submission by EDF .................................................................................................. 23
3.3
Submission by the Finnish Meteorological Institute (FMI) ..................................... 25
3.4
Submission by Idaho National Laboratory ............................................................... 26
INL Group 1 .......................................................................................................................... 26
INL Group 2 .......................................................................................................................... 27
3.5
Submission by IRSN ................................................................................................ 28
3.6
Submission by KAERI ............................................................................................. 29

4.

Observations ........................................................................................................................ 33
4.1
4.2

5.

Case 1 observations .................................................................................................. 33
Case 2 observations .................................................................................................. 34

Conclusions .......................................................................................................................... 37

References ................................................................................................................................... 39
Annex A. Submission by EDF ................................................................................................... 40
Executive summary .................................................................................................................. 40
A1.1 Introduction ...................................................................................................................... 45
A1.1.1. The OECD context ................................................................................................... 45
A1.1.2. Data .......................................................................................................................... 45
Notations ............................................................................................................................... 47
A1.2 Case 1 - Know model producing the synthetic data......................................................... 48
A1.2.1. The 3-points interpolation ........................................................................................ 48
A1.2.2. L2-error .................................................................................................................... 49
A1.2.3. L∞-error .................................................................................................................... 52
A1.2.4. Models comparison .................................................................................................. 54
A1.3 Case 2 - Unknown model producing the synthetic data................................................... 55
A1.3.1. Case 2a ..................................................................................................................... 55
A1.3.2. Case 2b..................................................................................................................... 60
A1.3.3. Case 2c ..................................................................................................................... 64
A1.4. Conclusion ...................................................................................................................... 68

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 7

Annex B. Submission by Finnish Meteorological Institute ..................................................... 69
Annex C. Submission by Idaho National Laboratory Group 1.............................................. 97
Introduction............................................................................................................................... 97
Assumptions ............................................................................................................................. 98
Assessment approach ................................................................................................................ 99
Results of the assessment........................................................................................................ 100
Annex D. Submission by Idaho National Laboratory Group 2............................................ 101
Introduction............................................................................................................................. 101
Assumptions ........................................................................................................................... 101
Model adequacy assessment approach.................................................................................... 101
Case 1 .................................................................................................................................. 102
Model adequacy assessment approach ................................................................................ 102
Results..................................................................................................................................... 106
Results of Assessment ......................................................................................................... 106
Case 2a ................................................................................................................................ 106
Assumptions ........................................................................................................................... 107
Model adequacy assessment approach.................................................................................... 107
Results..................................................................................................................................... 109
Results of assessment .......................................................................................................... 109
Case 2b ................................................................................................................................ 110
Results..................................................................................................................................... 112
Results of assessment .......................................................................................................... 113
Case 2c ................................................................................................................................ 113
Results..................................................................................................................................... 116
Results of assessment .......................................................................................................... 116
Annex E. Submission by IRSN ................................................................................................ 117
E.1.1. Case 1 – Known model producing the synthetic data .................................................. 117
1.
AM/GEV frequency model: ................................................................................... 118
2.
The POT GPD/Exponential frequency model: ....................................................... 120
3.
Case 1 – Conclusion ............................................................................................... 122
E1.2. Case 2 – Unknown model producing the synthetic data ............................................... 123
Cases 2a and 2b – Conclusion ............................................................................................. 124
Case 2a: the MA/GEV frequency model: ........................................................................... 124
Case 2b: the MA/GEV frequency model: ........................................................................... 126
Annex F. Submission by KAERI............................................................................................. 129
Introduction............................................................................................................................. 129
Synthetic data analysis ............................................................................................................ 129
F.1.1. Case 1 result assessment ........................................................................................... 131
F.1.2. Study for Case 2 ........................................................................................................ 132
F.1.3. Case 2 result assessment ........................................................................................... 135
Conclusion .............................................................................................................................. 136

List of figures
Figure 1.1. Example of the variety in potential modelling choices for magnitude-probability representation of a
streamflow initiating event
Figure 2.1. Case 1 plot up to a return period of 10 000 years
Figure 2.2. Case 2 plot up to a return period of 10 000 years
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

16
20
22

8 | NEA/CSNI/R(2021)10
Figure 3.1. EDF submission results plotted for the four cases
24
Figure 3.2. Fitting of the SM-Case1 data sets (w=100 years) with GPD distribution
28
Figure 3.3. To the left: Fitting Case2-a synthetic data with a GEV distribution; to the right: Fitting Case2a synthetic
data with a GEV distribution
29
Figure 3.4. KAERI submission results plotted for Case 1
31
Figure 3.5. KAERI submission results plotted for Case 2
32
Figure 4.1. Comparison of all submission results for Case 1
33
Figure 4.2. Comparison of all submission results for Case 2 for the 50 000 year predictions
35
Figure A1.1. Frequency - Magnitude obtained with the optimal GEV distribution (in blue, red and green) that best
predicts the known model or the synthetic data (in black) with respect to several error criteria
41
Figure A.1.2. Case 1: Frequency - Magnitude relation obtained with the optimal GEV distribution (in green) thatbest predicts the known model g (in black) with 3 points: in natural scale (left) and in logscale (right). Zoom on theinterval [1, 100] years on the first line.
50
Figure A1.3. Case 1: Relative error between the optimal GEV distribution and the given model g.
50
Figure A1.4. Case 1 Frequency - Magnitude relation obtained with the optimal GEV distribution (in red) that best
predicts the known model g (in black) with respect to the L2-error: in natural scale (left) and in logscale (right).
Zoom on the interval [1, 100] years on the first lines.
51
Figure A1.5. Case 1: Relative L2-error between the optimal GEV distribution and the given model g.
52
Figure A1.6. Case 1: Frequency - Magnitude relation obtained with the optimal GEV distribution (in red) thatbest predicts the known model g (in black) with respect to the L∞-error: in natural scale (left) and in logscale(right). Zoom on the interval [1, 100] years on the first line.
53
Figure A1.7. Case 1: Relative L∞-error between the optimal GEV distribution and the given model g
53
Figure A1.8. Case 1: Frequency - Magnitude obtained with the optimal GEV distribution (in blue, red and green)that best predicts the known model g (in black) with respect to the L2-error (in blue), the L∞-error (in red) and with 3points interpolation (in green): in natural scale (left) and in logscale (right).
54
Figure A1.9. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) with respect to the L2-error: in natural scale (left) and in logscale (right). 55
Figure A1.10. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in red) that best
predicts the given synthetic data (in black) with respect to the weighted L2-error: in natural scale (left) and in
logscale (right)
57
Figure A1.11. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) with respect to the L∞-error: in natural scale (left) and in logscale(right).
58
Figure A1.12. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) with respect to the weighted L∞-error: in natural scale (left) and in
logscale (right).
59
Figure A1.13. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue and red)that best predicts the given synthetic data (in black): in natural scale (left) and in logscale (right).
59
Figure A1.14. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) with respect to the L2-error: in natural scale (left) and in logscale(right).
60
Figure A1.15. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) thatbest predicts the given synthetic data (in black) with respect to the weighted L2-error: in natural scale (left)
and in logscale (right).
61
Figure A1.16. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) thatbest predicts the given synthetic data (in black) with respect to the L∞-error: in natural scale (left) and in logscale(right).
62
Figure A1.17. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) thatbest predicts the given synthetic data (in black) with respect to the weighted L∞-error: in natural scale (left) andin logscale (right).
63
Figure A1.18. Frequency - Magnitude relation obtained with the optimal GEV model (in blue and red) that best
predicts the given synthetic data (in black): in natural scale (left) and in logscale (right).
64
Figure A1.19. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) minimising the random mean squared error E Err2(Sc, qθ): in naturalscale (left) and in logscale (right).
65
Figure A1.20. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue) that best
predicts the given synthetic data (in black) minimising the Quantile0.95 Err2(Sc, qθ): in natural scale (left) andin logscale (right).
66

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 9
Figure A1.21. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue and red)that best predicts the given synthetic data (in black) minimising the mean feature of Err2(Sc, qθ) (in blue) and its
quantile 95% feature (in red): in natural scale (left) and in logscale (right).
Figure C.1. Generalised extreme value distribution family
Figure D.1. Return period versus mean magnitude
Figure D.2. Residuals vs fitted
Figure D.3. Normal Q-Q
Figure D.4. Scale-location
Figure D.5. Residuals vs leverage
Figure D.6. Magnitude vs return period
Figure D.7. Residuals
Figure D.8. Autocorrelation
Figure D.9. Normal Q-Q plot of standardised residuals
Figure D.10. Series redis(mod) [1:26]
Figure D.11. Series redis(mod) [1:26]
Figure D.12. Return period vs mean magnitude
Figure D.13. Residuals versus fitted
Figure D.14. Normal Q-Q
Figure D.15. Autocorrelation
Figure D.16. Principal component analysis results
Figure D.17. Return period versus mean magnitude
Figure D.18. Residuals versus fitted
Figure D.19. Normal Q-Q
Figure E.1. Fitting with a GEV distribution with the SM-Case1 data sets
Figure E.2. Fitting of the SM-Case1 data sets (w=100 years) with GPD distribution (u = 0.5 & λ = 1) for theare exponential).
return levels (the exceedances overFigure E.6. Hazard curve (Case 2b synthetic data)
Figure E.7. Fitting the Case 2b synthetic data
Figure E.8. Fitting Case 2a synthetic data with a GEV distribution
Figure F.1. Log regression fitting
Figure F.2. Fitting result for Case 1
Figure F.3. Magnitude prediction fitting for Case 1
Figure F.4. Linear regression fitting
Figure F.5. Linear regression fitting according to confidence interval
Figure F.6. Fitting result for Case 2
Figure F.7. Magnitude prediction fitting for Case 2

u

67
99
103
105
105
105
106
107
108
108
109
110
110
111
112
112
112
114
115
115
116
119
122
127
128
128
130
131
132
133
134
135
136

List of tables
Table 2.1. Synthetic data for Case 1
Table 2.2. Synthetic data for Case 2(a)Table 2.3. Synthetic data for Case 2(b)Table 2.4. Synthetic data for Case 2(c)Table 2.5. Results for the two cases
Table 3.1. Results for EDF submission (example of L2 criteria results)
Table 3.2. Case 1 results for FMI submission
Table 3.3. Case 2a results for FMI submission
Table 3.4. Case 2b results for FMI submission
Table 3.5. Case 2c results for FMI submission
Table 3.6. OpenBUGS script for Case 1 from INL Group 1 submission
Table 3.7. Case 1 results for INL Group 1 submission
Table 3.8. Case 2a results for INL Group 1 submission
Table 3.9. Case 2b results for INL Group 1 submission
Table 3.10. Case 2c results for INL Group 1 submission
Table 3.11. Case 1 results for INL Group 2 submission
Table 3.12. Case 2a results for INL Group 2 submission
Table 3.13. Case 2b results for INL Group 2 submission

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

20
21
21
21
22
24
25
25
25
25
26
26
26
26
27
27
27
27

10 | NEA/CSNI/R(2021)10
Table 3.14. Case 2c results for INL Group 2 submission
27
Table 3.15. Case 1 results for IRSN submission
29
Table 3.16. Case 2a results for IRSN submission
29
Table 3.17. Case 2b results for IRSN submission
29
Table 3.18. Case 1 results for KAERI submission
30
Table 3.19. Case 2 results for KAERI submission
30
Table 4.1. Case 1 results for all group submissions
34
Table 4.2. Case 2 results for all group submissions
34
Table A1.1. Case 1: Comparison between the known model g and the optimal GEV model obtained using several
error criteria
42
Table A1.2. Case 2a: Comparison between the given synthetic data and the optimal GEV model obtained using
several error criteria
42
Table A1.3. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, using several error criteria
43
Table A1.3. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, using several error criteria (Continued)
44
Table A1.4. Case 2c: Comparison between the given synthetic data and the optimal GEV model minimising
several features of the random error
44
Table A1.5. Case 1: synthetic data
45
Table A1.6. Case 2a: synthetic data
46
Table A1.7. Case 2b: additional synthetic data with respect to Case 2a46
Table A1.8. Case 2c: known return levels and uncertain ones
47
Table A1.9. Case 1: Comparison between the known model g and the optimal GEV model obtained with a
3-points interpolation (in %)49
Table A1.10. Case 1: Case 1: Comparison between the known model g and the optimal GEV model obtainedwith the relative L2-error (in %)52
Table A1.11. Case 1: Comparison between the known model g and the optimal GEV model obtained with the
L∞-relative error (in %)54
Table A1.12. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theL2-error, with the relative error in %.
56
Table A1.13. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theweighted L2-error, with the relative error (%).
57
Table A1.14. Case 2a: Comparison between the given synthetic data and the optimal GEV model using the
L∞-error, with the relative error (%).
58
Table A1.15. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theweighted L∞-error, with the relative error (%).
59
Table A1.16. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, with the relative error (%).
60
Table A1.17. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theweighted L2-error, with the relative error (%).
62
Table A1.18. Case 2b: Comparison between the given synthetic data and the optimal GEV model using the
L∞-error, with the relative error (%).
63
Table A1.19. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theweighted L∞-error, with the relative error (%).
64
Table A1.20. Case 2c: Comparison between the given synthetic data and the optimal GEV model minimising
the random mean squared error E Err2(Sc, qθ), with the relative error (%). For the uncertain values, we givethe quantile of the estimated return level.
65
Table A1.21. Case 2c: Comparison between the given synthetic data and the optimal GEV model minimising
the Quantile0.95 (Err2(Sc, qθ), with the relative error (%). For the uncertain values, we give the quantile of theestimated return level.
66
Table C.1. Synthetic data for Case 1.
97
Table C.2. Synthetic data for Case 2a98
Table C.3. Synthetic data for Case 2b
98
Table C.4. Synthetic data for Case 2c98
Table C.5. OpenBUGS script for Case 1 using GEV model
100
Table C.6. Case 1 results
100
Table C.7. Case 2a results
100
Table C.8. Case 2b results
100
Table C.9. Case 2c results
100
Table D.1. Resulting magnitude prediction for Case 1
106
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 11
Table D.2. Resulting magnitude prediction for Case 2-1
Table D.3. Resulting magnitude prediction for Case 2b
Table D.4. Benchmark data provided for Case 2-3
Table D.5. Resulting magnitude prediction for Case 2-3
Table E.1. Comparison of magnitudes (m) calculated with the synthetic data for Case 1. The values in brackets
correspond to the 70% confidence intervals.
Table E.2. The T-year quantiles with absolute and relative widths of their 70% confidence intervals
Table E.3. Synthetic data for Case 2a
Table E.4. Synthetic data for Case 2b
Table F.1. Synthetic data for Case 1
Table F.2. Regression results for Case 1
Table F.3. Optimisation for parameters
Table F.4. Comparison of magnitudes
Table F.5. Magnitude prediction by return period (Case 1)
Table F.6. Synthetic data for Case 2
Table F.7. Regression results for Case 2
Table F.8. Optimisation for parameters
Table F.9. Compare for magnitude
Table F.10. Magnitude prediction by return period (Case 2)

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

109
112
113
116
122
123
124
127
129
130
131
131
132
132
133
134
134
135

12 | NEA/CSNI/R(2021)10

List of abbreviations and acronyms
AM

Annual maxima

ANVS

Autoriteit Nucleaire Veiligheid en Stralingsbescherming (Netherlands)

CSNI

Committee on the Safety of Nuclear Installations (NEA)

EDF

Électricité de France

FMI

Finnish Meteorological Institute

GEV

Generalised extreme value

GPD

Generalised Pareto distribution

IE

Initiating event

INL

Idaho National Laboratory

IRSN

Institut de Radioprotection et de Sûreté Nucléaire (Institute forRadiological Protection and Nuclear Safety, France)

KAERI

Korea Atomic Energy Research Institute

MCMC

Markov Chain Monte Carlo

NEA

Nuclear Energy Agency

OECD

Organisation for Economic Co-operation and Development

POT

Peaks-Over-Threshold

PRA

Probabilistic risk assessment

SSE

Sum square error

TGNEV

Task Group on Natural External Events (NEA)

TNC

Truncated Newton Constrained

TSO

Technical and scientific support organisation

WGEVWorking Group on External Events (NEA)WGRISK

Working Group on Risk Assessment (NEA)

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 13

Executive summary
The March 2011 accident at the Fukushima Daiichi Nuclear Power Plant triggered
discussions about natural external events that are low-frequency but high-consequence.
To address these issues and determine which events would benefit from international cooperative work, the NEA Committee on the Safety of Nuclear Installations (CSNI)
established the Task Group on Natural External Events (TGNEV) at its June 2013
meeting. In June 2014, the CSNI decided to re-organise TGNEV into a Working Groupon External Events (WGEV) to improve the understanding and treatment of external
hazards and support the continued safety performance of nuclear installations as well as
improve the effectiveness of regulatory practices in NEA member countries. The WGEV
is composed of a forum of experts for the exchange of information and experience on
external events in NEA member countries, thereby promoting co-operation andmaintenance of an effective and efficient network of experts.
At its 61st meeting, the CSNI approved the recommended task on the benchmark on
external events hazard frequency and magnitude statistical modelling, to be pursued by
the WGEV. Modelling of these external events is a common practice in hazard and risk
assessments in many countries. Having a valid statistical approach to model thesehazards is important. However, current practice indicates a wide variety of approaches
being used and a lack of appreciation of the uncertainties inherent to these types ofstatistical models. The objective of this activity is to provide a benchmark suitable toexplore the application of typical approaches to external hazard representation througha data-informed process. This report captured two types of benchmarks, one with data
and model provided and one with just data provided (a “blind test”).
In the report is a summary of statistical modelling approaches from the organisations
Électricité de France (EDF), the Finnish Meteorological Institute (FMI), the Idaho
National Laboratory (INL), the Institut de Radioprotection et de Sûreté Nucléaire
(IRSN), and the Korea Atomic Energy Research Institute (KAERI).
Several observations can be made related to the approaches used and results from thebenchmark for external hazards.
•

Different statistical approaches such as regression or probability distribution
models can provide reasonable hazard frequency and magnitude estimations fortime periods where data exist.

•

Predictions for long return periods (e.g. much greater than the existing data timecollection) can prove challenging for some types of models and data sets.

•

Rather than focusing on predictions of a magnitude for a particular hazard, analternative approach might be to evaluate the probability of exceeding a criticallevel in a future time interval.

•

If an underlying physical phenomenon that drives an external hazard isunbounded then predictions for the hazard may be underestimated.

•

Capturing uncertainties in the hazard predictions is not typically performed.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

14 | NEA/CSNI/R(2021)10
Based upon the observations made from the benchmark, several potential future WGEV
activities are recommended:•

The availability of statistical open-source tools and frameworks offers thepotential for standardised approaches for representing the frequency/magnitude
hazards and should be investigated for hazard and risk applications.

•

Having knowledge of the underlying hazard phenomena could improve thepredictions made from models. For example, availability of a maximum upperbound (e.g. a physical limit) could help to improve long-term predictions. Future
benchmarks should investigate the application of physical phenomena withhazard modelling.

•

The use of “paleo-data” for a period longer than that recorded in actual data sets
could have resulted in better predictions from the statistical models. The
availability of this type of data, though, is not well understood for some types ofhazards, and should be investigated for external hazards of interest to the WGEV.

•

Uncertainties inherent in hazard model predictions should be better understoodand quantified as a part of validation.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 15

1. Introduction
1.1

BackgroundAs an input into risk analysis modelling and simulation, a hazard initiating event (IE) is
typically considered as the starting point for risk models. Since the IE both contributes
to the risk quantification results and provides the boundary conditions for the rest of thehazard scenarios, effective modelling of the frequency and magnitude for different
external events using data-driven methods can be applied. However, current practice has
shown a variety of technical approaches, models, and limitations in validation of theseapproaches (for example, see Figure 1.1). Consequently, this benchmark study isintended to demonstrate and capture commendable practices in formulating andassessing the quantification of external event IEs when using statistical models.
This benchmark was open to a variety of technical communities including academia,
government agencies, industry, research institutes, and technical and scientific support
organisations (TSOs). Chapter 2 of this report provides the case studies underconsideration of this benchmark. An overview of the results of the participant
submissions is presented in Chapter 3. Overall insights and conclusions are provided in
Chapter 4. Detailed submissions are listed in the Annexes.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

16 | NEA/CSNI/R(2021)10
Figure 1.1. Example of the variety in potential modelling choices for magnitude-probability representationof a streamflow initiating event

Source: OECD (2014), OECD Economic
http://dx.doi.org/10.1787/data-00688-en.

1.2

Outlook:

Statistics

andProjections(database),

Objective
The objective of this benchmark study was to facilitate an exercise on statisticalmodelling in order to better understand the quantitative technical analysis steps andprocesses used for assessing hazard frequency and magnitude in external events risk
assessments. This benchmark study report provides details (data and overall objectives)
for the benchmarking exercise (in Chapter 2) by specifying synthetic data for ahypothetical external event (e.g. precipitation, extreme temperatures, high winds).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 17

1.3

Overview of the exercises
In terms of scope, the benchmark exercise on the analysis and assessment of the hazard
frequency/magnitude (called “the benchmark” in this report) for external events risk
assessment covers the description of probabilistic hazard modelling, its associateduncertainty characterisation, and the process used for analysis and assessment. The
information that was provided to, and asked of, the benchmark participants
encompassed:•

Hypothetical observational data representing an external hazard that have beencreated from synthetic models (this type of model is used to create synthetic data
generated by a computer). Two cases are described: (1) a fully revealed “open”
case where both the synthetic data and the synthetic model producing the data
are provided, and (2) a “blind-test case” where only the synthetic data areprovided.

•

Descriptions by the participants of the assumptions made to create the hazard
frequency/magnitude model(s), the qualitative and quantitative results of themodel(s), the process used to assess the adequacy of the model(s) and the resultsof the model adequacy assessment. For this benchmark, phenomenologicallybased evaluation and modelling will not be considered since only “observationaltype” data (derived from a synthetic model) are provided, resulting in statisticaltypes of models to be considered. However, more complicated modelling
situations incorporating phenomena physics and/or spatial considerations may
be proposed for a future benchmark.

The specific choice of model(s) to be considered by benchmark participants is not limiteda priori but will need to be able to incorporate the data and should, ideally, be able toprovide a consideration of prediction uncertainty.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

18 | NEA/CSNI/R(2021)10

1.4

Glossary
This section contains a list of definitions for key terms used in this benchmark.
Aleatory uncertainty – Pertaining to stochastic (non-deterministic) events, the outcomeof which is described by a probability. From the Latin alea (game of chance, die).
Adequacy assessment – The process of judging, both qualitatively and quantitatively,
the predictive ability of a model used for decision making.
Data – Distinct observed (e.g. measured) values of a physical process. Generally, data
may be subject to uncertainties, such as imprecision in measurement.
External event – An event originating outside a nuclear power plant that directly or
indirectly causes an initiating event.
Hazard – Anything that has the potential to cause an undesired event or condition that
leads to damage.
Hazard curve – A model that relates the occurrence frequency of a hazard to themagnitude (e.g. intensity of an earthquake, precipitation rate, flood water level,temperature level) of the hazard.
Model – A mathematical construct that converts data and information into knowledge.
Two types of models are used for risk analysis purposes, probabilistic (or aleatory) anddeterministic.
Statistical model – A model that represents complex phenomena stochastically.
Examples of common statistical approaches used in risk assessment include extreme
value, exponential, Weibull and Poisson models.
Synthetic data – Data that are generated from computational models instead of from
actual observation.
Synthetic model – A computational model that is designed to produce data that mimic
actual observed data.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 19

2. Benchmark problem description
2.1

Data generation via synthetic models
Rather than real data gathered from the external hazard community, synthetic examples
of a hypothetical physical process (e.g. local precipitation, high winds) are used for thebenchmark. The reason for this is that because synthetic cases are used, the underlying
data-producing mechanisms are known (exactly), as opposed to real cases, where they
are not known. Thus, it is possible to assess how well data-driven models perform froma predictive standpoint. Furthermore, since synthetic models are created and run on acomputer to simulate data, it is possible to control the data-generating process, including
accounting for elements such as the uncertainty present in the synthetic data.
In general, a synthetic model is an equation that has the functional form of:Model output = f(inputs) * Uncertainty
For example, a simple linear model with no uncertainty used to represent a hypotheticalhazard frequency/magnitude relationship could be expressed as:
Magnitude = a * t
where the magnitude could be represented in terms of some observable quantity (e.g.
height about a river flood stage, quantity of rain in an hour, velocity of wind, depth ofsnow in a day), a is a changeable parameter used to control the output of the synthetic
model, and t is a time interval – or return frequency – that “sees” the observable quantity
described by the magnitude term. In this case, longer return frequencies would produce
larger synthetic events in a linear fashion.

2.2

Case 1 – Known model producing the synthetic data
The synthetic model used for the first exercise (Case 1) is:
M = 0.5 + 0.5 * log10(a * t).
This synthetic example is constructed so that different values of “return time intervals”
t produce a hypothetical (but known since it comes from the synthetic model) magnitude
M for an annual maxima event. It is possible to use these types of models to produce“synthetic data” where different event outcomes are predicted as a function of time (e.g.
producing a flooding hazard curve). As used in this benchmark, the synthetic model
serves as surrogate for a complex phenomenological process.
One characteristic in Case 1 is that the data generation process produces “extreme
events” such that the larger the event, the less frequently the event is seen – in otherwords the return interval becomes large as the event magnitude increases. Consequently,
the units (italicised) that are present in the Case 1 synthetic model are:M metres = 0.5 metres + 0.5*log10(a years-1 * t years) metres.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

20 | NEA/CSNI/R(2021)10
For Case 1, the “a variable” is set to 1. Thus, it is now possible to plot the magnitude ofthe event that is “produced” by Case 1 for select times up to a return interval of 10 000
years (see Figure 2.1).
Figure 2.1. Case 1 plot up to a return period of 10 000 years

However, what is important for this benchmark is that hypothetical observational data(from a synthetic model) is provided for the hazard frequency/magnitude modelling. For
Case 1, this data is shown in Table 2.1.
Table 2.1. Synthetic data for Case 1
Return period (years)

1

2

5

10

50

100

500

1 000

2 000

10 000Magnitude (metres)

0.50

0.65

0.85

1.0

1.4

1.5

1.9

2.0

2.2

2.5

Participants were asked to use the data for Case 1 and provide a model that best described
the frequency/magnitude relationship and the associated analysis and insights. The
results of this analysis should include those areas identified in Chapter 3 of thisbenchmark, including:•

Qualitative aspects and insights
o assumptions made to create the hazard frequency/magnitude model;
o the process used to assess the adequacy of the model.

•

Quantitative aspects and insights includingo the type of model describing the hazard frequency and magnitude statisticalresults;o uncertainties of the model - assessing uncertainty is important for both
validation and prediction (NRC, 2010);
o results of the model adequacy assessment or validation.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 21

2.3

Case 2 – Unknown model producing the synthetic data
For Case 2, the synthetic model used was not provided to the participants. However,three parts were provided for this example: part (a), which provides the synthetic data(ten data points) with no uncertainty on the data points provided; part (b), which provides
additional synthetic data (26 data points) with no uncertainty on the data points provided;and part (c), which has uncertainty estimates on some of the data. Any of the parts could
have been evaluated by participants.
For part (a), the synthetic data output from the unknown model are shown in Table 2.2.
Table 2.2. Synthetic data for Case 2(a)Return period (years)
Magnitude (metres)

1
0.53

2
0.53

5
0.54

10
0.55

50
0.59

100
0.62

500
0.79

1 000
0.95

3 000
1.6

10 0004.0

For part (b), the synthetic data output from the unknown model are shown in Table 2.3.
Table 2.3. Synthetic data for Case 2(b)Return period (years)
Magnitude (metres)
Return period (years)
Magnitude (metres)
Return period (years)
Magnitude (metres)

1
0.53
50
0.59
200
0.67

2
0.53
60
0.60
300
0.71

5
0.54
70
0.60
400
0.75

10
0.55
80
0.61
500
0.79

15
0.56
90
0.62
750
0.87

20
0.56
100
0.62
1 000
0.95

25
0.57
125
0.63
3 000
1.57

30
0.57
150
0.65
1 0000
3.97

40
0.58
175
0.66

For part (c), the synthetic data for long time intervals are presumed to not be known
exactly. For these times (500 years and longer), an estimate of the uncertainty on themagnitude has been provided and is shown in Table 2.4.
Table 2.4. Synthetic data for Case 2(c)Return period (years)
Magnitude (metres)

Mean
SDev.*
5th***
95th

1
0.53
-**
-

2
0.53
-

5
0.54
-

10
0.55
-

50
0.59

100
0.62
-

500
0.79
0.04
0.72
0.85

1 000
0.95
0.06
0.85
1.1

3 000
1.6
0.15
1.3
1.8

Note:
* SDev. is the standard deviation.
** A “-” indicates information that is not available.
*** The 5th and 95th indicates a 5th percentile and a 95th percentile value, respectively. In other words, there
is a 0.05 probability that the magnitude is less than or equal to 0.72 metres for an event with a return periodof 500 years.

The synthetic model used for Case 2 was:M metres = 0.5 metres + 0.01 * b * exp (1 years-1 * t0.19 years) metres.
where b is normal distribution with a mean of 1.1 and a standard deviation of 0.15. The
use of a normal distribution in this synthetic model provides a degree of “noise” orstochastic variation in the data points that are provided via the model. Using this model,
it is possible to plot the magnitude of the event that is “produced” by Case 2 for selecttimes up to a return interval of 10 000 years (see Figure 2.2).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

10 000
4.0
0.46
3.2
4.7

22 | NEA/CSNI/R(2021)10
Figure 2.2. Case 2 plot up to a return period of 10 000 years

Lastly, the magnitude of the model prediction was asked of the participants for 500,
5 000, 50 000, and 500 000 years.
Now that the two models are known completely, it is possible to provide the exact results
for the different return periods. These results are shown in Table 2.5. Note that since theCase 2 model uses an exponential, for long return periods, the magnitude can becomelarge – this model may not represent physical processes effectively but can provide achallenge for modelling.
Table 2.5. Results for the two cases
Return period (years)
Magnitude (metres) Case 1
Magnitude (metres) Case 2

500
1.9
0.78

5 000
2.4
2.2

50 0002.9
28

500 000
3.4
2 000

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 23

3. Benchmark outcomes
3.1

Benchmark completion
To complete the benchmark exercise, the following four areas should be addressed:
1. Assumptions
2. Model detail results
3. Model adequacy assessment approach
4. Results of the assessment
While much of the effort spent in responding to this benchmark was on technicalanalysis, it is important to understand key assumptions behind the modelling approachand the ultimate use of the model for hazard representation. As part of an integrated risk
analysis, the hazard characterisation is a key element of effective modelling. Thus, any
supporting analyses need to be developed in context with consideration andunderstanding of the scope, limitations, boundary conditions, complexity and contexts
that provide the foundation for the resulting statistical model. These types ofassumptions need to be clearly documented to understand how results will be used inrisk-informed decision making.
Submissions were provided for consideration by the following: Électricité de France(EDF), the Finnish Meteorological Institute (FMI), the Idaho National Laboratory (INL),
Institut de Radioprotection et de Sûreté Nucléaire (IRSN), and the Korea Atomic Energy
Research Institute (KAERI).

3.2

Submission by EDF
In its submission, EDF models the frequency-magnitude relation with a generalised
extreme value (GEV) distribution. It notes that, under some general hypotheses, theannual extreme value of a process, once normalised, tends to a GEV distribution and thereturn levels are specific quantiles of the annual extreme value distribution. The GEV
limit model makes it possible to estimate large return levels.
In Cases 1 and 2, the fitting is based on the minimisation of a criteria, defined as thesquared error, the weighted square error, the maximum error or the weighted maximum
error between the data and the model. For Case 1, EDF also used three-points
interpolation by solving a system of equations defined by the quantiles of GEV
distribution.
The EDF analysts noted it is possible to improve the precision of the GEV model onlarge return periods by penalising the errors to give more importance in properly
predicting the magnitude associated to large return periods.
In Case 2c, the data are uncertain. EDF uses a GEV model that minimises random
features of the error focusing on either the mean or 95th quantile. Examples of analysis
results obtained with the L2 criteria are shown in Table 3.1. Instead, the models were

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

24 | NEA/CSNI/R(2021)10
used to compare against the data tables that were shown in Chapter 2. Plots of theseresults are shown in Figure 3.1. The details of their analysis, including tables ofnumerical results, are shown in Annex A.
Table 3.1. Results for EDF submission (example of L2 criteria results)
Return period (years)
Case 1Case 2aCase 2bCase 2c

500
1.85
0.78
0.79
0.79

5 000
2.35
2.40
2.41
2.41

50 000
2.85
14.6
14.3
14.25

500 000
3.35
105
102
101

Figure 3.1. EDF submission results plotted for the four cases

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 25

3.3

Submission by the Finnish Meteorological Institute (FMI)The FMI submitted a write-up and included the analysis files (in the R tool). All thenecessary steps to repeat the exercise (excluding setting up some parts of the Renvironment) were provided in this notebook. The submission provides the theoreticaland technical background for the modelling approach and the assumptions. Then, theresults for fitting the selected statistical models are presented for each test case together
with an analysis of potential uncertainties.
The model used by FMI is the GEV. They use a Bayesian approach (using a couple ofnumerical approaches) to the estimation of the GEV parameters. One of the mainstrengths in the Bayesian approach is that it provides a natural way to estimate bothparameter and observational uncertainties. The numerical results of the FMI analysis areshown in 3.1 for Case 1 and 3.2 for Case 2.
Table 3.2. Case 1 results for FMI submission
Return period (years)
Magnitude (metres) Case 1 Exact
FMI mean from the Approximate Bayesian Computation
approach magnitude (metres)
FMI mean from traditional Markov Chain Monte Carlo approach(metres)

500
1.9
1.85

5 000
2.4
2.35

50 0002.9
2.84

500 000
3.4
3.34

1.84

2.35

2.88

3.43

Table 3.3. Case 2a results for FMI submission
Return period (years)
Magnitude (metres) Case 2 Exact
FMI mean from the Approximate Bayesian Computation
approach magnitude (metres)
FMI mean from traditional Markov Chain Monte Carlo approach(metres)

500
0.78
0.76

5 000
2.2
2.33

50 00028
16.62

500 000
2000
148.36

0.76

2.31

17.25

162.27

Table 3.4. Case 2b results for FMI submission
Return period (years)
Magnitude (metres) Case 2 Exact
FMI mean from the Approximate Bayesian Computation
approach magnitude (metres)
FMI mean from traditional Markov Chain Monte Carlo approach(metres)

500
0.78
0.77

5 000
2.2
2.32

50 00028
16.22

500 000
2000
142.11

0.76

2.30

16.90

156.36

Table 3.5. Case 2c results for FMI submission
Return period (years)
Magnitude (metres) Case 2 Exact
FMI mean from the Approximate Bayesian Computation
approach magnitude (metres)
FMI mean from traditional Markov Chain Monte Carlo approach(metres)

500
0.78
0.78

5 000
2.2
2.28

50 00028
14.16

500 000
2000
110.94

0.77

2.36

16.95

151.60BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

26 | NEA/CSNI/R(2021)10

3.4

Submission by Idaho National Laboratory
Two independent groups at the INL provided results. The first group used a Bayesian
approach for a GEV model, using the OpenBUGS analysis tool. The second group useda regression approach. Specifically, the second group used linear regression withtransformed form model to fit the magnitude vs. return period relationship model for
Case 1 and non-linear regression for Case 2.

INL Group 1
Group 1 used a Bayesian approach with the GEV model to quantify both Case 1 and 2.
The tool used was OpenBUGS. An example of the script used for Case 1 is shown in
Table 3.6. The predicted results for Case 1 are shown in Table 3.7. The results for Case
2 are shown in Table 3.8 (Case 2a), Table 3.9 (Case 2b), and Table 3.10 (Case 2c).
Table 3.6. OpenBUGS script for Case 1 from INL Group 1 submission

model{ for(i in 1:N) {z.p[i] ~ dnorm(mean[i],prec)
y.p[i] <- -log(1 - p[i])mean[i]<- mu - sigma/xi*(1 -pow(y.p[i],-xi))}
mu ~ dnorm(0,0.0001)
prec<-pow(sd,-2)
sd~dunif(0,10)xi ~ dunif(-1,1)
sigma ~ dunif(0,10)}
datalist(p=c(0.632, 0.393, 0.181, 0.0952, 0.0198, 0.00995, 0.002, 0.001, 0.0005, 0.0002, 0.0001, 0.00002,
0.000002),
z.p=c(0.50, 0.65, 0.85, 1.0, 1.4, 1.5, 1.9, 2.0, 2.2, NA, 2.5, NA, NA), N=13)
list(mu=1.0, sigma=1.0, xi=1.0)
Table 3.7. Case 1 results for INL Group 1 submission
Return period (years)
Magnitude (metres) Case 1 Exact
INL mean (metres)

500
1.9
1.88

5 000
2.4
2.37

50 0002.9
2.84

500 000
3.4
3.31

Table 3.8. Case 2a results for INL Group 1 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76

5 000
2.2
2.34

50 00028
16.02

500 000
2 000
136.20

Table 3.9. Case 2b results for INL Group 1 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76

5 000
2.2
2.32

50 00028
16.26

500 000
2 000
142.80BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 27

Table 3.10. Case 2c results for INL Group 1 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76

5 000
2.2
2.34

50 00028
16.02

500 000
2 000
136.20INL Group 2
The hazard frequency/magnitude model here for Case 1 fit the linear regression modelwith log transformation of the return period to describe the relationship between
magnitude and return period since the synthetic model provided is linear. The differenceis the residual term in regression equation, which is the vector values of the differences
between observed values and predicted values:� 𝚤𝚤 = 𝛽𝛽�0 + 𝛽𝛽�1 𝑙𝑙𝑙𝑙𝑙𝑙(𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝1 ) + 𝜖𝜖𝑖𝑖
𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀

The estimated equation is:� 𝚤𝚤 = 0.503364 + 0.506251 ∗ 𝑙𝑙𝑙𝑙𝑙𝑙(𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝1 )𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀

The estimated values above are the least square estimates of the intercept and slope. They
have standard error for the intercept and slope are 0.07006 and 0.03083, respectively.
The predicted results are shown in Table 3.11 for Case 1, Table 3.12 for Case 2a, Table
3.12 for Case 2b, and Table 3.14 for Case 2c.
Table 3.11. Case 1 results for INL Group 2 submission

Return period (years)
Magnitude (metres) Case 1 Exact
INL mean (metres)

500
1.9
1.87

5 000
2.4
2.38

50 0002.9
2.88

500 000
3.4
3.39

Table 3.12. Case 2a results for INL Group 2 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.75

5 000
2.2
2.32

50 00028
14.83

500 000
2 000
37.18

Table 3.13. Case 2b results for INL Group 2 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.75

5 000
2.2
2.28

50 00028
17.57

500 000
2 000
170.46

Table 3.14. Case 2c results for INL Group 2 submission
Return period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.71

5 000
2.2
2.01

50 00028
13.75

500 000
2 000
37.09

Group 2 noted that principal component regression has a few technical limitations anduncertainties. First, the uncertainty of mean for short time interval was not provided,even though the unknown standard deviation of each mean magnitude is small. To
estimate the standard deviation, the 5th and 95th percentile of the mean magnitude, an
imputation must be utilised which produced some uncertainty on those imputed values.
Compared with the results estimation obtained from Case 2a, we noticed that the
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

28 | NEA/CSNI/R(2021)10
confidence and prediction intervals for Case 2c are all wider than those for Case 2a.
Model fitting also does not perform well when there is a long time return period so thatthe predicted values do not fall into the 95% confidence interval. Lastly, the principal
component regression model could not meet the expected condition that the meanmagnitude keeps increasing as the return period increases.

3.5

Submission by IRSN
The IRSN submitted its analyses for Case 1, Case 2(a) and Case 2(b) in a write-up that
included the analysis files (using the R language).
For Case 1, a set of 100 return periods was first randomly sampled. Then, associatedmagnitudes were calculated using the synthetic model. Finally, extreme frequency
estimations were performed.
The analyst implemented an annual maxima (AM) frequency model in which the
distribution of the extreme events converges to a GEV one. The IRSN noticed that the
GEV distribution can have the form of the synthetic model only if the shape parameterof the distribution is equal to -1 and the scale parameter is negative. With suchconditions, the theoretical upper tail can only be finite and bounded and the GEV
distribution cannot be used with a negative scale parameter: thus it was not used todescribe the frequency/magnitude relationship for Case 1.
The analyst implemented a Peaks-Over-Threshold (POT) frequency model in which the
distribution of the exceedances over the threshold converge to an exponential one(Generalised Pareto Distribution, or GPD). The choice of the threshold value is based onusing the synthetic model. The Renext R library (developed by the IRSN and Alpestat)
was used for the frequency estimations. The relative difference in estimated magnitudes
does not exceed 5%, and all the plotting positions are inside the one sigma confidence
interval even though the latter is very narrow (see Figure 3.2).
Figure 3.2. Fitting of the SM-Case1 data sets (w=100 years) with GPD distribution

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 29

Table 3.15. Case 1 results for IRSN submission
Return period (years)
Magnitude (metres) Case 1 Exact
GPD/Exp IRSN mean (metres)

500
1.9
1.78

5 000
2.4
2.25

50 0002.9
2.72

500 000
3.4
3.19

For both Cases 2a and 2b, non-linear least-squares estimates of the GEV parameters
were performed. As shown in Figure 3.3, the fitting with confidence intervals is quite

good with heavy tails (very high shape parameter ξ ≈ 0.96 ) up to some thousands ofyears. For longer return periods the fitting is off the curve, perhaps because the proposedcases ignore any physical limits.
Figure 3.3. To the left: Fitting Case2-a synthetic data with a GEV distribution; to the right: Fitting
Case2a synthetic data with a GEV distribution

Table 3.16. Case 2a results for IRSN submission
Return period (years)
Magnitude (metres) Case 2 Exact
GEV IRSN mean (metres)

500
0.78
0.75

5 000
2.2
2.32

50 00028
16.79

500 000
2 000
150.34

Table 3.17. Case 2b results for IRSN submission
Return period (years)
Magnitude (metres) Case 2 Exact
GEV IRSN mean (metres)

3.6

500
0.78
0.76

5 000
2.2
2.31

50 00028
16.50

500 000
2 000
146.53

Submission by KAERIAs a result of fitting the relationship between magnitude and return period in Case 1 andCase 2, the relationship was estimated in log and linear regression. The data fitting ishigh in R2 but with a square error. The reason is sensitive to the coefficients of thesynthetic model. Therefore, further parameter analysis is required. The regression
equation of this study is as shown in the following equations (1) to (2).
For Case 1: 𝑀𝑀 = 𝐴𝐴 ∗ 𝑙𝑙𝑙𝑙 (𝑥𝑥) + 𝐵𝐵 (𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵)

For Case 2: 𝑀𝑀 = 𝐴𝐴 ∗ (𝑥𝑥) + 𝐵𝐵 (𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵)

(1)
(2)

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

30 | NEA/CSNI/R(2021)10
To minimise the square error of parameters A and B, an optimisation method is applied
through the square error solver function. As a result, the estimated parameters A and B
for Case 1 and Case 2 are shown in the following equations (3) and (4).
For Case 1: 𝑀𝑀 = 0.219861834 ∗ 𝑙𝑙𝑙𝑙 (𝑥𝑥) + 0.503363549 (𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒)(3)

For Case 2: 𝑀𝑀 = 0.000344252 ∗ (𝑥𝑥) + 0.565051348 (𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀 𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠𝑠 𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒) (4)The predicted results for Case1 and Case2 are shown in Table 3.17 and 3.18. In addition,
the re-fitted graph was compared with the existing data and displayed, as shown in Figure
3.4 to Figure 3.5.
Table 3.18. Case 1 results for KAERI submission
Return period (years)
Magnitude (metres) Case 1 Exact
Log KAERI mean (metres)
Optimisation KAERI mean (metres)

500
1.9
1.870
1.870

5 000
2.4
2.376
2.376

50 000
2.9
2.883
2.882

500 000
3.4
3.389
3.388

50 00028
15.56517.778500 000
2 000
150.565
172.691

Table 3.19. Case 2 results for KAERI submission
Return period (years)
Magnitude (metres) Case 2 Exact
Linear KAERI mean (metres)
Optimisation KAERI mean (metres)

500
0.78
0.715
0.737

5 000
2.2
2.065
2.286

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 31

Figure 3.4. KAERI submission results plotted for Case 1

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

32 | NEA/CSNI/R(2021)10
Figure 3.5. KAERI submission results plotted for Case 2

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 33

4. Observations
This benchmark provided an exercise focused on the quantitative technical analysis steps
and processes used to assess hazard frequency and magnitude for external events risk
assessments. The research provided details (data and overall objectives) for two uniquebenchmarking exercises specific to external events hazard frequency and magnitude
modelling by providing synthetic data for a hypothetical external event (precipitation,
extreme temperatures, high winds, etc.) to be determined by the participants. The
information and results provided by benchmark participants were collected andsummarised to gain insights on best practices from the activity.

4.1

Case 1 observations
Figure 4.1 shows the overall results of all participants for Case 1. As can be seen in the
figure, the results tended to track well with the exact results and were consistent acrossall groups. The numerical values used for Figure 4.1 are listed in Table 4.1.
Figure 4.1. Comparison of all submission results for Case 1

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

34 | NEA/CSNI/R(2021)10
Table 4.1. Case 1 results for all group submissions
Result
Exact
EDF
FMI
INL1
INL2
IRSN
KAERI

1
0.50

2
0.65
0.58

5
0.85
0.83

10
1.00
0.99

50
1.35
1.35

100
1.50
1.50

0.49
0.50
0.50
0.50

0.65
0.66

0.86
0.86

1.37
1.36

0.66

0.86

1.02
1.01
0.90
1.01

1.53
1.52
1.42
1.52

1.36

500
1.85
1.85
1.84
1.88
1.87
1.78
1.87

1 000
2.00
2.00

2 000
2.15
2.15

2.03
2.02

2.17
2.17

2.02

2.18

5 000
2.35
2.35
2.35
2.37
2.38
2.25
2.38

10 000
2.50
2.50
2.51
2.53
2.53

50 000
2.85
2.85
2.88
2.84
2.88
2.72
2.88

500 000
3.35
3.35
3.43
3.31
3.39
3.19
3.39

Note: EDF = Électricité de France (example of results obtained with the L2 criteria), FMI = Finnish
Meteorological Institute, INL1 = Idaho National Laboratory Group 1, INL2 = Idaho National Laboratory
Group 2, IRSN = Institut de Radioprotection et de Sûreté Nucléaire, KAERI = Korea Atomic Energy
Research Institute.
The approach and models used by the different groups ranged from curve fitting logarithm models toperforming Markov Chain Monte Carlo (MCMC) calculations using the GEV model. For Case 1, five of thegroups used a regression-based approach while two used a MCMC calculation.

4.2

Case 2 observations
The study for Case 2 was much more complicated since the underlying model was notprovided. A summary of the numerical predictions that were provided by the submitters
is shown in Table 4.2. As can be seen in the table, the predictions for long periods oftime (> 10 000 years) started to become problematic due to the underlying synthetic
model (based upon an exponential function that was proposed as a strong challenge).
Table 4.2. Case 2 results for all group submissions
Result
EDF
FMI
INL1
INL2
IRSN

Exact

Case
2a
2b
2c
2a
2b
2c
2a
2b
2c
2a
2b
2c
2a
2b
2c

500
0.78
0.78
0.79
0.79
0.76
0.76
0.77
0.76
0.76
0.76
0.75
0.75
0.71
0.75
0.76

Return period (years)
5 000
50 000
2.20
28.0
2.40
14.6
2.41
14.3
2.41
14.25
2.31
17.3
2.30
16.9
2.36
17.0
2.34
16.0
2.32
16.3
2.34
16.0
2.32
14.8
2.28
17.6
2.01
13.8
2.32
16.8
2.31
16.5

500 000
2 000
105
102
101
162
156
152
136
143
136
37
170
37
150
147

Note: EDF example of results obtained with the L2 criteria

Focusing on just the 50 000-year magnitude predictions (see Figure 4.2), it is possible tosee that the predictions were consistent, but too low by about 40%. Again, this under
prediction is likely a result of the use of the exponential function for the underlying
synthetic model for Case 2. While the models used by the submitters (e.g. GEV,Generalised Pareto, regression) can represent exponential types of behaviour, the data
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 35

that was provided may represent too short of a period to reliably predict very long return
periods (i.e. for periods of time much longer than the data set provided). In Case 2, the“longest” data point provided was for 10 000 years, or five time shorter than the
50 000-year prediction.
Case 2b represented the situation where additional data was provided (over that provided
by Case 2a) – 26 data points were provided instead of the original ten. However, in bothcases, the longest time period represented was 10 000 years. Looking at the 50 000-year
predictions, we can see that increasing the data points (from 10 to 26) resulted in only
slightly better predictions (but they were still low by approximately 40%). It appears thathaving additional data within the seen time did not lead to much of an improvement in
the predictions. Alternatively, if we had provided “paleo-data” for a period longer than
10 000 years, this data would have likely resulted in better predictions. This aspect of
the analysis might be considered for future benchmark exercises.
Figure 4.2. Comparison of all submission results for Case 2 for the 50 000 year predictions

The participants that submitted information for Case 2 used two basic approaches to theanalysis, Bayesian analysis via MCMC or regression to a curve. It appears that thequantity of data and the quantification had little impact on the results shown in
Figure 4.2. However, it did appear that the Bayesian approach was able to quantify theuncertainty on the predictions in a straight-forward fashion (for example, the estimate
for the INL1 2c result was a mean of 16 m with a 5th percentile of 14 m and a 95th
percentile of 17 m for the 50 000-year return period).
When we look at the 500 000-year magnitude predictions, we can see that the predictions
were consistent, but very low. This under prediction is (again) likely a result of the useof the exponential function for the underlying synthetic model for Case 2. Note that theuse of an exponential function was chosen as a challenging case and most likely does
not represent any real physical phenomena.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

36 | NEA/CSNI/R(2021)10
However, a concern here would be if return periods much greater than what has been
experienced (> 10 times the time period) were being predicted and the underlying
phenomena that drive the magnitude of the external hazard is unbounded upward like an
exponential function. In that case, the predictions for the hazard may be underestimated.
Fortunately, most (if not all) natural external hazards would have upper bounds to themagnitude – it would not make sense to assume infinitely large floods, temperatures,rain, snow, earthquakes, etc.
The difficulty in the predictions for long return periods does indicate that knowledge ofthe underlying phenomena or knowledge of upper bounds (e.g. physical limits) could
help to greatly improve the predictions. This aspect of the analysis might be considered
for future benchmark exercises.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 37

5. Conclusions
This benchmark study provided a mechanism to develop statistics-based approaches forrepresenting complex physical phenomena such as extreme temperatures, droughts, highwinds, floods and extreme snowfall. The benchmark provided the overall objectives ofthe study, the data to be investigated and the expectations for results reporting.
This report provided two independent benchmark cases to facilitate understanding ofstatistical modelling that may be used in quantitative technical analysis to assess hazard
frequency and magnitude for external events risk assessment:•

Case 1: Based on a data set with a given associated synthetic model (i.e. known
function form).

•

Case 2: Based on a data set generated from another associated synthetic modelthat is not given (i.e. a “blind study”).

Participants were asked to provide any key assumptions, the statistical model(s) used,
the overall results of the analysis, and the process used to assess the adequacy of themodel(s), including uncertainties. The provided data were created using synthetic
models for a hypothetical external event (e.g. precipitation, extreme temperatures andhigh winds). As such, the data used in this benchmark may not be representable for any
real physical phenomenon and are solely intended to test different modelling approaches.
One of the findings of this benchmark study is that different statistical approaches, suchas regression or probability distribution model application, can provide reasonablehazard frequency and magnitude estimation for time periods where data exist (i.e.
interpolation) for both cases. However, the benchmark study did identify potential
issues:•

Predictions for long return periods (much greater than the existing data timecollection) can prove challenging for some types of data sets.

•

Rather than focusing on predicting a magnitude for a particular hazard, analternative approach might be to evaluate the probability of exceeding a criticallevel in a future time interval.

•

If an underlying physical phenomenon that drives an external hazard isunbounded (e.g. like the exponential function used in Case 2), predictions for thehazard may be underestimated.

•

Some participants had a process to quantify uncertainties. However, capturing
the uncertainties in the predictions was not typically performed.

While issues and challenges do exist in statistical modelling of hazards, this study
pointed to potential improvements in hazard modelling processes. It is recommended
that the following items be considered for future WGEV activities:

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

38 | NEA/CSNI/R(2021)10•

The availability of statistical open-source tools and frameworks offers thepotential for standardised approaches for representing the frequency/magnitudeof hazards applicable to risk applications.

•

While not part of this benchmark, having knowledge of the underlying hazard
phenomena could improve the predictions made from models. For example, theavailability of a maximum upper bound (e.g. a physical limit) could help togreatly improve long-term predictions.

•

The use of “paleo-data” for a period longer than that recorded in actual data sets
could have resulted in better predictions from the statistical models. The
availability of this type of data, though, is not well understood for some types ofhazards.

•

Uncertainties inherent in hazard model predictions should be better understoodand quantified as a part of validation.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 39

References
National Research Council (2010), Assessing the Reliability of Complex Models,
https://www.nap.edu/read/13395/chapter/1
OECD (2014), OECD Economic Outlook: Statistics and Projections (database),
http://dx.doi.org/10.1787/data-00688-en.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

40 | NEA/CSNI/R(2021)10

Annex A. Submission by EDF
Executive summary
The modelling of external hazards encompasses different technical aspects, depending on
the type of hazard. All these hazards impact nuclear installations. Usually, a specific risk
analysis is performed for a given hazard to inform a decision. A risk analysis contains a setof scenarios, frequencies and associated consequences. A scenario is a sequence that
contains an initiating event and one or more subsequent events. The end state (aconsequence) of a scenario is also of interest.
The global objective of the benchmark study launched by the OECD NEA, is to focus on
an initiating event induced by an external hazard. The statistical modelling of a hazard
frequency and magnitude is of particular interest in this benchmark.
The benchmark study provides synthetic data for a hypothetical external event.
Hypothetical data describing the external hazard comes from synthetic models. These data
have been generated with a computer code. Two cases are described:•

Case 1: a fully revealed open case where both the synthetic data and the synthetic
model producing the data are provided;•Case 2: a blind-test case where only the synthetic data are provided.
In this report, we model the frequency-magnitude relation with a GEV distribution.
Indeed, under some general hypotheses, the annual extreme value of a process, once
normalised, tends to a GEV distribution. Furthermore, the return levels are specificquantiles of the annual extreme value distribution. The GEV limit model makes it
possible to estimate large return levels.
To fit the GEV model that best predicts the given synthetic model or the given
synthetic data, we used several criteria:•

the minimum of the squared error between the data and the model;•

the maximum error over a large range of the return period.
In addition, for Case 1, we used a three-points interpolation by solving a system of
equations defined by the quantiles of GEV distribution. We voluntarily present thismethod only in Case 1 but it can be applied on the other cases without uncertainty.
It is possible to improve the precision of the GEV model on large return periods. To
this end, we penalise the errors to give more importance to properly predicting themagnitude associated with large return periods.
If the data are uncertain, we define a GEV model that minimises random features ofthe error: the mean or its quantile of order 95%.
For each case, we present the relative error and the optimal set of parameters of the
GEV model.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 41

Figure A1.1 draws the T -return level obtained with the optimal GEV distribution in eachcase.
Table A1.1 to Table A1.4 detail the reference magnitudes and the approximated ones,
obtained using several error criteria.
Figure A1.1. Frequency - Magnitude obtained with the optimal GEV distribution (in blue, red and green)that best predicts the known model or the synthetic data (in black) with respect to several error criteria

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

42 | NEA/CSNI/R(2021)10
Table A1.1. Case 1: Comparison between the known model g and the optimal GEV model obtained using
several error criteria
Returnperiod(years)Ref
Mag(m)ApproxMag(m) 3pinterpApprox
Mag(m) L2error
Approx
Mag(m) L8error

2

5

10

50

100

500

1 000

2 000

10 000

0.650

0.849

1.0

1.349

1.5

1.849

2.0

2.151

2.5

0.582

0.828

0.99

1.349

1.5

1.85

2.0

2.151

2.5

0.600

0.841

1.0

1.353

1.503

1.850

2.0

2.150

2.50

0.724

0.901

1.02

1.31

1.44

1.75

1.90

2.04

2.40

Table A1.2. Case 2a: Comparison between the given synthetic data and the optimal GEV model obtainedusing several error criteria
Returnperiod(years)Ref Mag(m)
ApproxMag (m) L2-error
ApproxMag (m) weighted
L2-error
ApproxMag (m) L8-error
ApproxMag (m) weighted
L8-error

2

5

10

50

100

500

1 000

3 000

10 000

0.53

0.54

0.55

0.59

0.62

0.79

0.95

1.6

4.0

0.523

0.526

0.530

0.556

0.584

0.776

0.986

1.730

3.959

0.505

0.507

0.511

0.537

0.566

0.759

0.971

1.719

3.969

0.479

0.482

0.486

0.511

0.540

0.732

0.942

1.685

3.915

0.445

0.448

0.452

0.478

0.507

0.701

0.914

1.669

3.940

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 43

Table A1.3. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, using several error criteria
ReturnPeriod(years)Ref Mag(m)
ApproxMag (m)- L2error
ApproxMag (m)weighted
L2-error
ApproxMag (m)- L8error
ApproxMag (m)weighted
L8-error
ReturnPeriod(years)Ref Mag(m)
ApproxMag (m)- L2error
ApproxMag (m)weighted
L2-error
ApproxMag (m)- L8error
ApproxMag (m)weighted
L8-error
ReturnPeriod(years)Ref Mag(m)
ApproxMag (m)- L2error

2

5

10

50

100

500

1 000

3 000

10 000

0.53

0.54

0.55

0.59

0.62

0.79

0.95

1.6

4.0

0.543

0.546

0.550

0.576

0.604

0.795

1.00

1.741

3.952

0.538

0.541

0.545

0.571

0.599

0.790

1.00

1.74

3.96

0.479

0.482

0.486

0.512

0.540

0.732

0.942

1.685

3.916

0.448

0.451

0.455

0.481

0.510

0.704

0.917

1.67

3.94

15

20

25

30

40

60

70

80

90

0.56

0.56

0.57

0.57

0.58

0.60

0.60

0.61

0.62

0.554

0.557

0.560

0.564

0.570

0.582

0.587

0.593

0.599

0.548

0.552

0.555

0.558

0.564

0.576

0.582

0.588

0.593

0.489

0.493

0.496

0.499

0.506

0.518

0.523

0.529

0.535

0.459

0.462

0.465

0.469

0.475

0.487

0.493

0.499

0.504

125

150

175

200

300

400

750

0.63

0.65

0.66

0.67

0.71

0.75

0.87

0.618

0.631

0.643

0.656

0.704

0.750

0.902

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

44 | NEA/CSNI/R(2021)10
Table A1.3. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, using several error criteria (Continued)
ReturnPeriod(years)ApproxMag (m)weighted
L2-error
ApproxMag (m)- L8error
ApproxMag (m)weighted
L8-error

2

5

10

50

100

500

1 000

0.612

0.626

0.638

0.651

0.699

0.745

0.897

0.554

0.567

0.580

0.592

0.641

0.687

0.840

0.523

0.537

0.550

0.562

0.611

0.658

0.813

3 000

10 000

Table A1.4. Case 2c: Comparison between the given synthetic data and the optimal GEV model
minimising several features of the random error
ReturnPeriod(years)Ref
Mag(mean)(m)
ApproxMag (m)- MeanApproxMag (m)Quantile
95%

2

5

10

50

100

500

1 000

3 000

10 000

0.53

0.54

0.55

0.59

0.62

0.79

0.95

1.6

4.0

0.525

0.527

0.531

0.558

0.587

0.781

0.993

1.73

3.94

0.517

0.519

0.523

0.551

0.581

0.778

0.994

1.75

3.98

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 45

A1.1 Introduction
A1.1.1. The OECD context
The modelling of external hazards encompasses different technical aspects depending on
the type of hazard, all of them impacting nuclear installations. There is also a link
between the hazard and the associated risk analysis. A risk analysis contains a set of
scenarios, frequencies and associated consequences, developed in such a way as toinform decisions. A scenario contains an initiating event and (usually) one or more
subsequent events leading to an end state that reflects the issue of concern.
The objective of the benchmark study launched by the OECD NEA is to focus on the IE
by facilitating an exercise on the statistical modelling for assessing hazard frequency andmagnitude for external event risk assessment. The benchmark study provides synthetic
data for a hypothetical external event.
The analysis steps and modelling results we obtained are detailed on this report.
We would like to underline that the scope of this benchmark differs from the activities
that EDF usually performs regarding extreme natural events characterisation. Indeed,
what is usually available is a set of measured data for a given phenomenon, and what is
performed is statistical extrapolation of these data to evaluate a magnitude of the natural
phenomenon for a high return level period.

A1.1.2. Data
Hypothetical observational data represent an external hazard which has been createdfrom synthetic models (this type of model is used to create synthetic data that have been
generated from a computer). Two cases are described:•

Case 1: a fully revealed open case where both the synthetic data and the synthetic
model producing the data are provided;•Case 2: a blind-test case where only the synthetic data are provided.

Case 1 - Known model producing the synthetic data
The synthetic model used for the first exercise is:g : T → g(T ) = 0.5 + 0.5 log10(T ) (1)where T is the return period in years and M = g(T ) the associated magnitude given in
metres. This model led to the return period / magnitude detailed in Table A1.5.
Table A1.5. Case 1: synthetic data
Return period (years)
Mag (m)

1

2

5

10

50

0.50

0.65

0.85

1.0

1.4

100

500

1 000

1.9

Case 2 - Unknown model producing the synthetic data
In this case, only the synthetic data are provided (the synthetic model used is notprovided). The Case 2 presents three subcategories:

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

46 | NEA/CSNI/R(2021)10•

On the first two, the synthetic data have no uncertainty and the last has
uncertainty estimates on some of the data.

•

The three-point interpolation would also be applicable for the 2a case data set(without uncertainty). In order to focus on presenting the differences between adata set with and without uncertainties, we focused the Case 2 studies on the L2and L∞-error models.

Case 2a: The synthetic data (ten data points) with no uncertainty on the points areprovided and gathered in Table A1.6.
Table A1.6. Case 2a: synthetic data
Return period (years)25

Mag (m)0.53

0.54 0.55

1050
10 000
0.59
4.0

100

500

1 000

3 000

0.62

0.79

0.95

1.6

Case 2b: In addition to the synthetic data of Case 2a, 16 data points are provided without
uncertainty. The additional data are seven.
Table A1.7. Case 2b: additional synthetic data with respect to Case 2a
Return period (years)
Mag

15
0.56

20
0.56

25
0.57

30
0.57

40
0.58

60
0.60

70
0.60

80
0.61

Return period (years)
Mag (m)

90
0.62

125
0.63

150
0.65

175
0.66

200
0.67

300
0.71

400
0.75

750
0.87

Care: In Table A1.7, we changed the value given in Table 2.3 of the benchmark(benchmark numeration). Case 2b: Comparison between the given synthetic data and the
optimal GEV model using the L2-error, using several error criteria of the benchmark(benchmark numeration) at the period T = 3 000 years from the value 1.57 to the value
1.6 and the value given in the benchmark at the period T = 10 000 years from the value
3.97 to the value 4.0, in order to make the data Table A1.7. Case 2b: Comparison
between the given synthetic data and the optimal GEV model using the L2-error, using
several error criteria (benchmark numeration) coherent with the data of Table A1.6. Case
2a: Comparison between the given synthetic data and the optimal GEV model obtainedusing several error criteria (benchmark numeration).
Case 2c: Here, the synthetic data for long time intervals (500 years and longer) are
presumed to not be known exactly. For these return periods, the uncertainty on themagnitude is provided. The other points (≤ 500 years) have no uncertainty on theassociated magnitude. The synthetic data (ten data points) are gathered in Table A1.8.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 47

Table A1.8. Case 2c: known return levels and uncertain ones
Return period (years)
Mean
St. dev.
q0.05
q0.95

1
0.53
-

2
0.53
-

5
0.54
-

10
0.55
-

50
0.59
-

100
0.62
-

500
0.79
0.04
0.72
0.85

1 000
0.95
0.06
0.85
1.1

3 000
1.6
0.15
1.3
1.8

10 000
4.0
0.46
3.2
4.7

We can note that the quantiles q0.05 and q0.95 given in Table A1.11 correspond to those of
a normal distribution for which the mean and standard deviation are those given in thetable.
Thus, we will model the uncertainty by a normal distribution as follows:
X500 ∼ Normal(µ = 0.79, σ = 0.04)(2)

X1000 ∼ Normal(µ = 0.95, σ = 0.06)(3)

X3000 ∼ Normal(µ = 1.6, σ = 0.15)(4)

X10 000 ∼ Normal(µ = 4.0, σ = 0.465)(5)

We consider that these four random variables are independent.

Notations
Let us note Gθ the cumulated density function of a GEV distribution para- metered by θ
= (µ, σ, ξ). The Extreme Value Theory shows that the annual maximum of theunderlying process follows a GEV distribution.
The T -return level qθ(T ) is by definition the quantile of order (1 − 1/T ) of the GEV
distribution :G θ(q θ(T )) = 1 − 1/T(6)

The T -return level is evaluated from the relation:(7)

In this report, we want to fit the GEV distribution that best predicts the return levels
given either by a known 1model (Case 1) or by some synthetic data (Case 2).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

48 | NEA/CSNI/R(2021)10

A1.2 Case 1 - Know model producing the synthetic data
In Case 1, the model g that produced the synthetic data is given by (1).
The objective of the study is to find the optimal parameter θ = (µ, σ, ξ) so that the GEV
distribution parametered by θ fits the best frequency-magnitude relation (1).
To find the optimal θ, we first considered a 3-points interpolation according to the relation(37). Then, we consider a criteria to evaluate the model error between g and qθ, noted
Err(g, qθ) and we find θ∗ defined by :(8)

In our study, we tested different criteria evaluating the model error Err(g, qθ):•the L2-error based criteria, which cumulates the errors at each return period T:
see Section A1.2.2;•

the L∞-error based criteria, which focuses on the maximal error over the
return period T : see Section A.1.2.3.

A1.2.1. The 3-points interpolation
The GEV distribution being defined by three parameters, we consider a system based on
three equations as follows:

where (Ti, Mi)i=1,2,3 represents the return period and its associated magnitude taken in
the Table A1.5.
To determine the parameter θ =( µ, σ, ξ), solution of this system, the parameters µ and σ
are defined as a function of ξ : µ = 𝑓𝑓 1 (𝑇𝑇 𝑀𝑀 ) (ξ) and σ = 𝑓𝑓 2 (𝑇𝑇 𝑀𝑀 ) (ξ), with 𝑓𝑓 1 and 𝑓𝑓 2
𝑖𝑖,

𝑖𝑖

𝑖𝑖,

𝑖𝑖

known. Then, we reach the parameter ξ, solution of the equation 𝑓𝑓 3 (𝑇𝑇 𝑀𝑀 ) (ξ) = 0, using𝑖𝑖,

𝑖𝑖

the Brent’s algorithm. Note that the functions 𝑓𝑓 1, 𝑓𝑓 2 and 𝑓𝑓 3 are presented in Annex A.
Once the optimal ξ∗ is obtained, we deduce µ∗ =𝑓𝑓 1 (𝑇𝑇 𝑀𝑀 ) (ξ*) and σ∗= 𝑓𝑓 2 (𝑇𝑇 𝑀𝑀 ) (ξ*).
𝑖𝑖,

𝑖𝑖

𝑖𝑖,

𝑖𝑖

The optimal parameter with (Ti, Mi)i=1,2,3 = {(10 000, 2.5) , (100, 1.5) , (2 000, 2.2)} is:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.502162, 0.216915, −1.1 · 10−6) (9)These points were chosen because they lead to a better approximation after several tests
performed in the Annex A.
With ξ∗ near to 0, we can say that the hazard frequency-magnitude relation follows a
Gumbel distribution with parameters (µ∗, σ∗) = (0.502162, 0.216915).
Figure A1.2 draws the T -return level obtained with the optimal GEV distribution andthe model g.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 49

Table A1.9 details the T -return levels obtained with the optimal GEV for several periods
T and gives the relative error made for each of them (in %) defined by:(10)

Figure A1.3 draws the relative error function: T  ε(T).
Table A1.9. Case 1: Comparison between the known model g and the optimal GEV model obtained with a
3-points interpolation (in %)Return period (years)
Ref Mag (m)Approx Mag (m)
ε(T )

2
0.650
0.582
10.58

5
0.849
0.828
2.57

10
1.0
0.99
9.7e-1

50
1.349
1.349
6.94e-2

100
1.5
1.5
0

500
1.849
185
2.68e-2

1 000
2.0
2.0
2.21e-2

2 000
2.151
2.151
1.55e-2

10 0002.5
2.5
0A1.2.2. L2-error
We define the error as the square of the L2-norm of (g − qθ). In other words, we consider
the integrated squared error at each point T between the model g(T ) and qθ(T ):(11)
To compute the integration, we use the Gauss-Legendre algorithm to approximate theintegral with a finite sum of the integrand evaluated on some judicious points. Using the
TNC (Truncated Newton Constrained) optimisation algorithm with 100 different starting
points, we obtain the following optimal point:θ∗ = (µ∗, σ∗ , ξ ∗ ) = (0.522748, 0.211424, 0.0032985)

The associated error is:(12)

Err∗ = 0.1356

As the 3-points interpolation, the hazard frequency-magnitude relation is defined
following the Gumbel distribution with parameters θ∗ = (µ∗, σ∗) = (0.522748,
0.211424).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

50 | NEA/CSNI/R(2021)10
Figure A.1.2. Case 1: Frequency - Magnitude relation obtained with the optimal GEV distribution (ingreen) that best predicts the known model g (in black) with 3 points: in natural scale (left) and in logscale(right). Zoom on the interval [1, 100] years on the first line.

Figure A1.3. Case 1: Relative error between the optimal GEV distribution and the given model g.

optimisation algorithm with 100 different starting points, we obtain the following
optimal point:θ∗ = (µ∗, σ∗ , ξ ∗ ) = (0.522748, 0.211424, 0.0032985)

The associated error is:(12)

Err∗ = 0.1356

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 51

As the 3-points interpolation, the hazard frequency-magnitude relation is defined
following the Gumbel distribution with parameters θ∗ = (µ∗, σ∗) = (0.522748,
0.211424).
Figure A1.4 draws the T -return level obtained with the optimal GEV distribution andthe model g.
(%)

Table A1.5 details the T -return levels obtained with the optimal GEV for several periods
T and the relative error made for each of them (in %).
Figure A1.6 draws the relative error function: T → ε(T ).
Figure A1.4. Case 1 Frequency - Magnitude relation obtained with the optimal GEV distribution (in red)that best predicts the known model g (in black) with respect to the L2-error: in natural scale (left) and in
logscale (right). Zoom on the interval [1, 100] years on the first lines.

Table A1.10. Case 1: Comparison between the known model g and the optimal GEV model obtained with therelative L2-error (in %)Return period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.650
0.600
7.72

5
0.849
0.841
1.04

10
1.0
1.0
3.0e-2

50
1.349
1.353
2.64e-1

100
1.5
1.503
1.83e-1

500
1.849
1.850
2.83e-2

1 000
2.0
2.0
6.43e-3

2 000
2.151
2.150
2.29e-2

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

10 0002.5
2.50
3.74e-3

52 | NEA/CSNI/R(2021)10
Figure A1.5. Case 1: Relative L2-error between the optimal GEV distribution and the given model g.

A1.2.3. L∞-error
We define the error as the L∞-norm of (g − qθ). In other words, we consider the
maximum absolute error between g(T ) and qθ(T ) :(13)

To solve the optimisation problem, we use the Cobyla algorithm initialised with the
optimal point obtained with the L2-error.
We obtain the following optimal point:θ∗ = (µ∗, σ∗, ξ∗ ) = (0.668508, 0.149587, 0.0484294)(14)

The associated error is:
Err∗ = 0.107711

Figure A1.6 draws the T -return level obtained with the optimal GEV distribution andthe model g.
Table details the T -return level obtained with the optimal GEV for several periods T
and gives the relative error made for each of them (in %) defined by (10).
Figure A1.7 draws the relative error function: T → ε(T ).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 53

Figure A1.6. Case 1: Frequency - Magnitude relation obtained with the optimal GEV distribution (in red)that best predicts the known model g (in black) with respect to the L∞-error: in natural scale (left) and in
logscale (right). Zoom on the interval [1, 100] years on the first line.

Figure A1.7. Case 1: Relative L∞-error between the optimal GEV distribution and the given model g

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

54 | NEA/CSNI/R(2021)10
Table A1.11. Case 1: Comparison between the known model g and the optimal GEV model obtained with
the L∞-relative error (in %)Return period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.650
0.724
11.3

5
0.849
0.901
6.09

10
1.0
1.02
2.41

50
1.349
1.31
2.85

100
1.5
1.44
4.05

500
1.849
1.75
5.22

1 000
2.0
1.90
5.22

2 000
2.151
2.04
5.0

10 000
2.5
2.40
3.81

A1.2.4. Models comparison
In order to ease comparison, Figure A1.8 draws both models obtained with 3-points
interpolation and with both criteria’s (11) and (13).
Figure A1.8. Case 1: Frequency - Magnitude obtained with the optimal GEV distribution (in blue, red andgreen) that best predicts the known model g (in black) with respect to the L2-error (in blue), the L∞-error(in red) and with 3-points interpolation (in green): in natural scale (left) and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 55

A1.3 Case 2 - Unknown model producing the synthetic data
For Case 2, the synthetic data is provided without the associated model.

A1.3.1. Case 2a
In Case 2a, we use the synthetic data of Table A1.6.

A1.3.1.1. L2-error
We define the following error criteria:(15)
where Sa = (Ti, xTi )1≤i≤10 are the synthetic data of Table A1.6, with Ti the return period
and xTi the associated magnitude.
Using the TNC optimisation algorithm with 100 different starting points, we got thefollowing optimal point:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.522859, 0.001, 0.868784)(16)

The associated error is:
Err∗ = 0.1517

Figure A1.9 draws the T -return level obtained with the optimal GEV distribution andthe sample Sa.
Figure A1.9. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the L2-error: in natural scale (left)
and in logscale (right).

Table details the T -return level obtained with the optimal GEV for several return
periods T and gives the relative error made for each of them (in %) defined by:

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

56 | NEA/CSNI/R(2021)10(17)

Table A1.12. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theL2-error, with the relative error in %.
Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.523
1.27

5
0.54
0.526
2.60

10
0.55
0.530
3.66

50
0.59
0.556
5.79

100
0.62
0.584
5.75

500
0.79
0.776
1.76

1 000
0.95
0.986
3.84

3 000
1.6
1.730
8.08

10 000
4.0
3.959
1.03

A1.3.1.2. Weighted L2-error
In order to force a best adequation over the large return periods, we consider theweighted least square error defined by:(18)2

where Sa = (Ti, xTi )1≤i≤10 are the synthetic data of Table A1.6, with Ti the return period
and xTi the associated magnitude.
Using this weighted error, we penalise the errors on large return periods, considering,for example, that we prefer getting a model that predicts better return levels associatedto large return periods even if the model works worse on low return periods.
Using the TNC optimisation algorithm with 100 different starting points, we obtain thefollowing optimal point:θ∗ = (µ∗, σ∗ , ξ ∗ ) = (0.504188, 0.001, 0.869813)(19)

The associated error is:
Err∗ = 0.4018

Figure A1.10 draws the T -return level obtained with the optimal GEV distribution andthe sample Sa.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 57

Figure A1.10. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in red)that best predicts the given synthetic data (in black) with respect to the weighted L2-error: in naturalscale (left) and in logscale (right)

Table A1.13 details the T -return level obtained with the optimal GEV for several periods
T and gives the relative error made for each of them (in %) defined by (17). We verify
that although the global error is larger, the optimal GEV distribution fits better the datafor periods larger than 500 years.
Table A1.13. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theweighted L2-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.505
4.79

5
0.54
0.507
6.06

10
0.55
0.511
7.06

50
0.59
0.537
8.94

100
0.62
0.566
8.73

500
0.79
0.759
3.95

1 000
0.95
0.971
2.17

3 000
1.6
1.719
7.44

10 000
4.0
3.969
0.78

A1.3.1.3. L∞-error
We define the following error criteria:
Err(Sa , qθ ) = max |xTi − qθ (Ti )|(20)

where Sa = (Ti, xTi )1≤i≤10 are the synthetic data of Table A1.6, with Ti the return period
and xTi the associated magnitude.
Using the TNC optimisation algorithm with the optimal point obtained with the L2-error
as a starting point, we obtain the following optimal point:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.478532, 0.00100012, 0.868792)

The associated error is :(21)

Err∗ = 0.0851144

Figure A1.11 draws the T -return level obtained from the optimal GEV distribution andthe sample Sa.
Table A1.14 details the T -return level obtained from the optimal GEV for several
periods T and gives the relative error made for each of them (in %).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

58 | NEA/CSNI/R(2021)10
Figure A1.11. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the L∞-error: in natural scale (left)
and in logscale (right).

Table A1.14. Case 2a: Comparison between the given synthetic data and the optimal GEV model using the
L∞-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.479
9.63

5
0.54
0.482
10.8

10
0.55
0.486
11.7

50
0.59
0.511
13.3

100
0.62
0.540
12.9

500
0.79
0.732
7.36

1 000
0.95
0.942
8.12e-1

3 000
1.6
1.685
5.32

10 000
4.0
3.915
2.12

A1.3.1.4. Weighted L∞-error
In order to force a best adequation over the large return periods T, we consider theweighted error defined byErr(Sa , qθ ) = max | log(Ti)[xTi − qθ (Ti )]|(22)

where Sa = (Ti, xTi )1≤i≤10 are the synthetic data of Table 1.2, with Ti the return period
and xTi the associated magnitude.
Using this weighted error, we penalise the errors on large T , considering, for example,
that we prefer getting a model that predicts better large return levels even if the model
works worse on little return levels.
Using the TNC optimisation algorithm with the optimal point obtained with the L2-error
as starting point, we obtain the following optimal point:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.44483, 0.001, 0.870939)

The associated error is:(23)

Err∗ = 0.553827

Figure A1.12 draws the T -return level obtained from the optimal GEV distribution andthe sample Sa.
Table A1.15 details the T -return level obtained from the optimal GEV for several
periods T and gives the relative error made for each of them (in %) defined by (17). We
note that weighting the error for large periods has changed very little the optimal

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 59

solution: it means that the maximum error was reached for large T , so that thepenalisation was not efficient.
Figure A1.12. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the weighted L∞-error: in naturalscale (left) and in logscale (right).

Table A1.15. Case 2a: Comparison between the given synthetic data and the optimal GEV model using theweighted L∞-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.445
15.98

5
0.54
0.448
17.05

10
0.55
0.452
17.85

50
0.59
0.478
18.98

100
0.62
0.507
18.26

500
0.79
0.701
11.28

1 000
0.95
0.914
3.764

3 000
1.6
1.669
4.318

10 000
4.0
3.940
1.479

A1.3.1.5. Models comparison
In order to ease comparison, Figure A1.13 draws both models obtained from both criteria(15), (18), (20) and (22).
Figure A1.13. Case 2a: Frequency - Magnitude relation obtained with the optimal GEV model (in blue
and red) that best predicts the given synthetic data (in black): in natural scale (left) and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

60 | NEA/CSNI/R(2021)10

A1.3.2. Case 2b
In Case 2b, we use the synthetic data of Table A1.6 and the additional ones of Table
A1.7.

A1.3.2.1. L2-error
We define the following error criteria:(24)2where Sb = (Ti , xTi )1≤i≤26 are the synthetic data of Tables A1.6 and A1.7, with Ti thereturn period and xTi the associated magnitude.
Using the TNC optimisation algorithm with 100 different starting points, we got thefollowing optimal point:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.542937, 0.001, 0.867794)

The associated error is:(25)

Err∗ = 0.1715

Figure A1.14 draws the T -return level obtained from the optimal GEV distribution andthe sample Sb.
Table A1.16 details the T -return level obtained from the optimal GEV for several
periods T and gives the relative error made for each of them (in %) defined by (17).
Figure A1.14. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the L2-error: in natural scale (left)
and in logscale (right).

Table A1.16. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theL2-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.543
2.52
15
0.56
0.554
1.16
125
0.63
0.618
1.97

5
0.54
0.546
1.11
20
0.56
0.557
5.44e-1
150
0.65
0.631
2.98

10
0.55
0.550
1.68e-2
25
0.57
0.560
1.71
175
0.66
0.643
2.51

50
0.59
0.576
2.40
30
0.57
0.564
1.14
200
0.67
0.656
2.10

100
0.62
0.604
2.55
40
0.58
0.570
1.76
300
0.71
0.704
8.19e-1

500
0.79
0.795
6.23e-1
60
0.60
0.582
3.04
400
0.75
0.750
4.17e-2

1 000
0.95
1.00
5.68
70
0.60
0.587
2.08
750
0.87
0.902
3.65

3 000
1.6
1.741
8.82
80
0.61
0.593
2.76

10 000
4.0
3.952
1.21
90
0.62
0.599
3.43

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 61

A1.3.2.2. Weighted L2-error
In order to force a best adequation over the large periods T, we consider the weighted
least square error defined by:

2

(26)

where Sb = (Ti, xTi )1≤i≤26 are the synthetic data of Table A1.6 and Table A1.7, with Ti
the return period and xTi the associated magnitude.
Using the TNC optimisation algorithm with 100 different starting points, we obtain thefollowing optimal point:θ∗ = (µ∗, σ∗, ξ ∗ ) = (0.53762, 0.001, 0.868286) (27)The associated error is:
Err∗ = 0.6837

Figure A1.15 draws the T -return level obtained from the optimal GEV distribution andthe sample Sb.
Table A1.17 details the T -return level obtained from the optimal GEV for several periods
T and gives the relative error made for each of them (in %) defined by (17). We note thateven if the global error is larger, the optimal GEV distribution fits better for periods
larger than 500 years.
Figure A1.15. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the weighted L2-error: in naturalscale (left) and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

62 | NEA/CSNI/R(2021)10
Table A1.17. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theweighted L2-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.538
1.52
15
0.56
0.548
2.11
125
0.63
0.612
2.79

5
0.54
0.541
1.30e-1
20
0.56
0.552
1.49
150
0.65
0.626
3.77

10
0.55
0.545
9.83e-1
25
0.57
0.555
2.63
175
0.66
0.638
3.29

50
0.59
0.571
329
30
0.57
0.558
2.07
200
0.67
0.651
2.86

100
0.62
0.599
3.39
40
0.58
0.564
2.67
300
0.71
0.699
1.52

500
0.79
0.790
2.99e-2
60
0.60
0.576
3.92
400
0.75
0.745
6.0e-2

1 000
0.95
1.00
5.25
70
0.60
0.582
2.96
750
0.87
0.897
3.15

3 000
1.6
1.74
8.74
80
0.61
0.588
3.62

10 000
4.0
3.96
1.00
90
0.62
0.593
4.27

A1.3.2.3. L∞-error
We define the error criteria defined in (20).
Err(Sb , qθ ) = max |xTi − qθ (Ti )|(28)

where Sb = (Ti, xTi )1≤i≤26 are the synthetic data of Table A1.6 and Table A1.7 with Ti
the return period and xTi the associated magnitude.
Using the TNC optimisation algorithm with the optimal point obtained with the L2-error
as starting point, we obtain the following optimal point:θ∗ = (µ∗, σ∗, ξ∗ ) = (0.478596, 0.00100039, 0.868782)(29)

The associated error is:
Err∗ = 0.0854126

Figure A1.16 draws the T -return level obtained from the optimal GEV distribution andthe sample Sb.
Table A1.18 details the T -return level obtained from the optimal GEV for several
periods T and gives the relative error made for each of them (in %) defined by (17).
Figure A1.16. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the L∞-error: in natural scale (left)
and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 63

Table A1.18. Case 2b: Comparison between the given synthetic data and the optimal GEV model using the
L∞-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.479
9.62
15
0.56
0.489
12.6
125
0.63
0.554
12.1

5
0.54
0.482
10.8
20
0.56
0.493
12.0
150
0.65
0.567
12.8

10
0.55
0.486
11.71
25
0.57
0.496
13.0
175
0.66
0.580
12.2

50
0.59
0.512
13.3
30
0.57
0.499
12.4
200
0.67
0.592
11.6

100
0.62
0.540
12.9
40
0.58
0.506
12.8
300
0.71
0.641
9.77

500
0.79
0.732
7.35
60
0.60
0.518
13.72
400
0.75
0.687
8.39

1 000
0.95
0.942
7.99e-1
70
0.60
0.523
12.8
750
0.87
0.840
3.50

3 000
1.6
1.685
5.34
80
0.61
0.529
13.3

10 000
4.0
3.916
2.10
90
0.62
0.535
13.8

A1.3.2.4. Weighted L∞-error
In order to force a best adequation over the large return periods T , we consider theweighted least square error defined by :Err(Sb , qθ ) = max | log(Ti)[xTi − qθ (Ti )]|(30)

where Sb = (Ti, xTi )1≤i≤26 are the synthetic data of Table A1.6 and Table A1.7, with Ti
the return period and xTi the associated magnitude.
Using the TNC optimisation algorithm with the optimal point obtained with the L2-error
as starting point, we got the following optimal point:θ∗ = (µ∗, σ∗ , ξ ∗ ) = (00.447999, 0.001, 0.870727)

The associated error is:(31)

Err∗ = 0.570137

Figure A1.17 draws the T -return level obtained from the optimal GEV distribution andthe sample Sb.
Table A1.19 details the T -return level obtained from the optimal GEV for several periods
T and gives the relative error made for each of them (in %) defined by (17). We note thateven if the global error is larger, the optimal GEV distribution fits better for periods
larger than 1 000 years.
Figure A1.17. Case 2b: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) with respect to the weighted L∞-error: in naturalscale (left) and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

64 | NEA/CSNI/R(2021)10
Table A1.19. Case 2b: Comparison between the given synthetic data and the optimal GEV model using theweighted L∞-error, with the relative error (%).
Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)e(T )Return Period (years)
Ref Mag (m)Approx Mag (m)
e(T )

2
0.53
0.448
15.4
15
0.56
0.459
18.1
125
0.63
0.523
16.9

5
0.54
0.451
16.5
20
0.56
0.462
17.5
150
0.65
0.537
17.4

10
0.55
0.455
17.3
25
0.57
0.465
18.3
175
0.66
0.550
16.7

50
0.59
0.481
18.4
30
0.57
0.469
17.8
200
0.67
0.562
16.1

100
0.62
0.510
17.8
40
0.58
0.475
18.1
300
0.71
0.611
13.9

500
0.79
0.704
10.9
60
0.60
0.487
18.8
400
0.75
0.658
12.2

1 000
0.95
0.917
3.49
70
0.60
0.493
17.8
750
0.87
0.813
6.59

3 000
1.6
1.67
4.41
80
0.61
0.499
18.2

1 0000
4.0
3.94
1.54
90
0.62
0.504
18.6

A1.3.2.5. Models comparison
In order to ease comparison, Figure A1.18 draws both models obtained from all thecriteria.
Figure A1.18. Frequency - Magnitude relation obtained with the optimal GEV model (in blue and red)that best predicts the given synthetic data (in black): in natural scale (left) and in logscale (right).

A1.3.3. Case 2c
In Case 2c, we use the synthetic data of Table A1.8.
To fit a GEV distribution on the synthetic data, we define the error function:2(32)

where Sc = (Ti, xTi )1≤i≤6 are the known synthetic data and (XTi )1≤i≤4 are the 4 unknown
return levels, only described by their distribution.
Then, the error function Err2(Sc, qθ) is a random variable because of the randomness of
the return levels associated to the return periods 500, 1 000, 3 000 and 10 000 years.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 65

We will consider several probabilistic features of the random variable Err2(Sc, qθ) tofind the optimal θ.

A1.3.3.1. Mean feature
In this section, we consider the mean of the random variable Err2(Sc, qθ) and we findthe optimal θ that minimises the mean error:θ∗ = argmin E Err2(Sc, q θ)(33)

θ with the optimal point obtained in the Case 2b
Using the TNC optimisation algorithmwith the L2-error as a starting point, we obtain the following optimal point:θ∗ = (µ∗, σ ∗, ξ∗) = (0.524142, 0.00104785, 0.862527)
The associated error is:(34)

Err∗ = 0.6969

Figure A1.19 draws the T -return level obtained with the optimal GEV distribution.
Table A1.20 details the T -return level obtained from the optimal GEV for all the periods
T. When the return level was unknown, we give to which quantile of the uncertainty
distribution the model value corresponds.
Figure A1.19. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) minimising the random mean squared error E
Err2(Sc, qθ): in natural scale (left) and in logscale (right).

Table A1.20. Case 2c: Comparison between the given synthetic data and the optimal GEV model
minimising the random mean squared error E Err2(Sc, qθ), with the relative error (%). For the uncertain
values, we give the quantile of the estimated return level.
Return Period (years)
Mean
Approx Mag (m)Quantile
e(%)

2
0.53
0.525
1.02

5
0.54
0.527
2.34

10
0.55
0.531
3.38

50
0.59
0.558
5.41

100
0.62
0.587
5.30

500
0.79
0.781
0.41
-

1 000
0.95
0.993
0.76
-

3 000
1.6
1.73
0.82
-

10 000
4.0
3.94
0.45
-

A1.3.3.2. Quantile feature
In this section, we consider the quantile 95% of the random variable Err2(Sc, qθ) and wefind the optimal θ that minimises the quantile error:
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

66 | NEA/CSNI/R(2021)10θ∗ = argmin Quantile0.95 (Err2(Sc, qθ))(35)

θUsing the truncated newton method
optimisation algorithm with the optimal point
obtained in the Case 2b with the L2-error as a starting point, we obtain the following
optimal point:θ∗ = (µ∗, σ∗ , ξ ∗ ) = (0.516115, 0.00108486, 0.859784) (36)The associated error is:
Err∗ = 1.04575
It means that P[ Err2(Sc, qθ) ≥ Err∗2 ] = 0.05.

Figure A1.20 draws the T -return level obtained from the optimal GEV distribution.
Table A1.21 details the T –return level obtained from the optimal GEV for all the periods
T.
When the return level is unknown, we give to which quantile of the uncertainty distribution the model
value corresponds.

Figure A1.20. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue)that best predicts the given synthetic data (in black) minimising the Quantile0.95 Err2(Sc, qθ): in naturalscale (left) and in logscale (right).

Table A1.21. Case 2c: Comparison between the given synthetic data and the optimal GEV model
minimising the Quantile0.95 (Err2(Sc, qθ), with the relative error (%). For the uncertain values, we givethe quantile of the estimated return level.
Return Period (years)
Mean
Approx Mag (m)Quantile
e(%)

2
0.53
0.517
2.53

5
0.54
0.519
3.81

10
0.55
0.523
4.80

50
0.59
0.551
6.61

100
0.62
0.581
6.33

500
0.79
0.778
0.39
-

1 000
0.95
0.994
0.77
-

3 000
1.6
1.75
0.84
-

10 000
4.0
3.98
0.48
-

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 67

A1.3.3.3. Models comparison
In order to ease comparison, Figure A1.21 draws both models obtained when minimising
the mean error and its quantile of order 95%. We note that both criteria lead to the same
optimal solution.
Figure A1.21. Case 2c: Frequency - Magnitude relation obtained with the optimal GEV model (in blue
and red) that best predicts the given synthetic data (in black) minimising the mean feature of Err2(Sc, qθ)(in blue) and its quantile 95% feature (in red): in natural scale (left) and in logscale (right).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

68 | NEA/CSNI/R(2021)10

A1.4. Conclusion
The benchmark study proposed an exercise on the statistical modelling. The main goal
is to better understand the quantitative technical analysis steps and processes used forassessing hazard frequency and magnitude for external event risk assessments. To thisend, we have the synthetic data for a hypothetical external event.
In this report, we modelled the hazard frequency-magnitude relation with a generalised
extreme value distribution. Several criteria to fit the GEV distribution: L2-error and L∞error. For the Case 1, a 3-points interpolation is performed.
For the synthetic data with no uncertainty, we saw that the criterion allowing the best
fitting of GEV distribution is L2-error. For the synthetic data with the uncertainty in somepoints, minimising the random mean squared error has the best precision.
In this benchmark, whatever the type of data, we managed to approximate the synthetic
model/data by a GEV distribution.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 69

Annex B. Submission by Finnish Meteorological Institute

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

70 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 71

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

72 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 73

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

74 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 75

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

76 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 77

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

78 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 79

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

80 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 81

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

82 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 83

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

84 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 85

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

86 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 87

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

88 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 89

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

90 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 91

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

92 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 93

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

94 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 95

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

96 | NEA/CSNI/R(2021)10

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 97

Annex C. Submission by Idaho National Laboratory Group 1
Introduction
This paper responds to an OECD Nuclear Energy Agency (NEA) report that calls for abenchmark exercise on the analysis and assessment of the hazard frequency andmagnitude for external events risk assessments. The objective of this benchmark study
is to focus on the statistical modelling of an external hazard initiating event (IE),
assessing its frequency and magnitude, that could be used in probabilistic risk
assessment (PRA) of external hazards. The NEA report specifies two cases ofhypothetical observational data created from synthetic models (which are used to create
synthetic data that have been generated from a computer) for a hypothetical external
event (e.g. precipitation, extreme temperatures, high winds): (1) a fully revealed “open”
case where both the synthetic data and the synthetic model producing the data areprovided, and (2) a “blind-test case” where only the synthetic data are provided. This
paper documents the author’s analysis steps, assumptions, insights and modelling results
for the exercise.
Case 1 – Known model producing the synthetic data
The following synthetic model is provided in Case 1 to create hypothetical observational
data for the hazard frequency/magnitude modelling:M = 0.5 + 0.5*log10(a*t)
The synthetic model serves as surrogate for a complex phenomenological process. It isof the form that different values of “return time intervals” t produce a hypothetical (butknown since it comes from the synthetic model) magnitude M for an annual maxima
event. These types of models can be used to produce “synthetic data” and predict
different event outcomes as a function of time (e.g. producing a flooding hazard curve).
The hypothetical observational data for the magnitude M (in metres) and return period t(in years) from the synthetic model with the “a variable” set to 1 is shown in Table C.1.
Table C.1. Synthetic data for Case 1.
Return Period (years)
Magnitude (metres)

1
0.50

2
0.65

5
0.85

10
1.0

50
1.4

100
1.5

500
1.9

1 000
2.0

2 000
2.2

10 0002.5

Participants are asked to use the data for Case 1 and provide a model that best described
the frequency/magnitude relationship and the associated analysis and insights. The
results of this analysis should include those areas identified in Chapter 3 of thisbenchmark, including:•oQualitative aspects and insights.
o Assumptions made to create the hazard frequency/magnitude model.
The process used to assess the adequacy of the model.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

98 | NEA/CSNI/R(2021)10•

Quantitative aspects and insights including:o The type/form of the model describing the hazard frequency and magnitude
statistical results.
o Uncertainties of the model. Assessing uncertainty is important for both
validation and prediction [NRC, 2010].
o Results of the model adequacy assessment or validation.

Case 2 – Unknown model producing the synthetic data
In Case 2, the synthetic model used is not provided. Instead, three sets of synthetic data
output from the unknown model are provided.
•Case 2a provides the synthetic data (ten data points) with no uncertainty
associated with the data points (see Table C.2).

•

Case 2b provides additional synthetic data (26 data points) with no uncertainty
associated with the data points (see Table C.3).

•

Case 2c has uncertainty estimates on some of the data (see Table C.4).
Table C.2. Synthetic data for Case 2a

Return Period (years)
Magnitude (metres)

1
0.53

2
0.53

5
0.54

10
0.55

50
0.59

100
0.62

500
0.79

1 000
0.95

3 000
1.6

10 000
4.0

20
0.56
100
0.62
1 000
0.95

25
0.57
125
0.63
3 000
1.57

30
0.57
150
0.65
10 000
3.97

40
0.58
175
0.66

1 000
0.95
0.06
0.85
1.1

3 000
1.6
0.15
1.3
1.8

10 000
4.0
0.46
3.2
4.7

Table C.3. Synthetic data for Case 2b
Return Period (years)
Magnitude (metres)
Return Period (years)
Magnitude (metres)
Return Period (years)
Magnitude (metres)

1
0.53
50
0.59
200
0.67

2
0.53
60
0.60
300
0.71

5
0.54
70
0.60
400
0.75

10
0.55
80
0.61
500
0.79

15
0.56
90
0.62
750
0.87

Table C.4. Synthetic data for Case 2c
Return Period (years)
Magnitude (metres)

Mean
Sdev.*
5th***
95th

1
0.53
**

2
0.53

5
0.54

10
0.55

50
0.59

100
0.62

500
0.79
0.04
0.72
0.85

Assumptions
A generalised extreme value (GEV) model was used in this benchmark exercise topredict the magnitude values for long return periods using the synthetic data provided inCase 1 and Case 2 (Table C.2., C.3. and C.4.). OpenBUGS was used to implement the
GEV model. It is believed by many analysts that observed data from relatively shortterm return periods should not be used to predict the magnitude values for much longerterm periods (e.g. use observed data up to 5 000 years to predict for 50 000 years).
However, when such needs arise, for example for the use of external hazard frequency
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 99

in probabilistic risk assessment (PRA) models, various approaches and engineering
judgement have to be used to provide such estimates, and the GEV model is believed tobe such a tool to provide best estimates, with the modelling uncertainty being in mind
when using the results from the model.

Assessment approach
The GEV is a single family of limiting distribution that combines three other limit
distributions: Gumbel, Fréchet and Weibull. It has the following cumulative distributionfunctional form:
1𝑧𝑧−𝜇𝜇 − �𝜉𝜉��}𝜎𝜎𝐺𝐺(𝑧𝑧) = exp {− �1 + 𝜉𝜉 �

(1)where 𝜉𝜉 is a shape parameter, 𝜇𝜇 is a location parameter, and 𝜎𝜎 is a scale parameter. The
shape parameter 𝜉𝜉 determines the distribution type and the tail behaviour. 𝜉𝜉 = 0
corresponds to the Gumbel distribution, which is unbounded and has exponential tail.
𝜉𝜉 > 0 corresponds to the Fréchet distribution, which has lower bound and long tail. 𝜉𝜉 <0 corresponds to the Weibull distribution, which has an upper bound and short tail(Figure C.1).
Figure C.1. Generalised extreme value distribution family

The inverse distribution function, or quantile function, is often used to calculate theextreme values:
𝜎𝜎𝑧𝑧𝑝𝑝 = 𝐺𝐺 −1 (1 − 𝑝𝑝) = 𝜇𝜇 − �1 − [− log(1 − 𝑝𝑝)]−𝜉𝜉 �, 𝜉𝜉 ≠ 0
𝜉𝜉(2)

The return level 𝑧𝑧𝑝𝑝 is exceeded by the annual maximum in any particular year with theannual exceedance probability 𝑝𝑝, when the measuring time is in years. Or one can saythat the return level 𝑧𝑧𝑝𝑝 is exceeded, on average, once during the return period 1/𝑝𝑝.
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

100 | NEA/CSNI/R(2021)10
The GEV model shown in Eq. (1) was applied to the synthetic data in this benchmark
exercise using OpenBUGS. An example of the script used for Case 1 is shown in Table
C.5.
Table C.5. OpenBUGS script for Case 1 using GEV model

model{ for(i in 1:N) {z.p[i] ~ dnorm(mean[i],prec)
y.p[i] <- -log(1 - p[i])mean[i]<- mu - sigma/xi*(1 -pow(y.p[i],-xi))}
mu ~ dnorm(0,0.0001)
prec<-pow(sd,-2)
sd~dunif(0,10)xi ~ dunif(-1,1)
sigma ~ dunif(0,10)}
datalist(p=c(0.632, 0.393, 0.181, 0.0952, 0.0198, 0.00995, 0.002, 0.001, 0.0005, 0.0002, 0.0001, 0.00002,
0.000002),
z.p=c(0.50, 0.65, 0.85, 1.0, 1.4, 1.5, 1.9, 2.0, 2.2, NA, 2.5, NA, NA), N=13)
list(mu=1.0, sigma=1.0, xi=1.0)

Results of the assessment
The predicted results for Case 1 are shown in Table C.6. The results for Case 2 are shownin Table C.7(Case 2a), Table C.8(Case 2b), and Table C.9(Case 2c).
Table C.6. Case 1 results
Return Period (years)
Magnitude (metres) Case 1 Exact
INL mean (metres)

500
1.9
1.88

5 000
2.4
2.37

50 0002.9
2.84

500 000
3.4
3.31

5 000
2.2
2.34

50 00028
16.02

500 000
2 000
136.20

5 000
2.2
2.32

50 00028
16.26

500 000
2 000
142.80

5 000
2.2
2.34

50 00028
16.02

500 000
2 000
136.20

Table C.7. Case 2a results
Return Period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76
Table C.8. Case 2b results

Return Period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76
Table C.9. Case 2c results

Return Period (years)
Magnitude (metres) Case 2 Exact
INL mean (metres)

500
0.78
0.76

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 101

Annex D. Submission by Idaho National Laboratory Group 2Introduction
The aim of the benchmark is to investigate the frequency/magnitude relationship forexternal events. Two cases represent different external hazard from synthetic data. The
first case provides both the synthetic data and synthetic model. The second caseintroduces three parts with only the synthetic data, where the first and second parts do
not consider the uncertainty, whereas the third part of the magnitude has been provided
with the uncertainty. The magnitude from those two cases is modelled as a function of
return period. The return period is also known as recurrence interval, which is theestimated average time between floods events occurring.
Linear regression with a transformed form model would be considered to fit themagnitude vs. return period relationship model for case1. Cases using ordinary leastsquared method would have the same assumptions and similar model adequacy
assessment approach. A linear model has the following general form:
�0 + 𝛽𝛽�1 𝑥𝑥1 + 𝛽𝛽
�2 𝑥𝑥2 … + 𝜖𝜖𝑖𝑖
𝑌𝑌𝑖𝑖 = 𝛽𝛽where the Y is the response variable, x are predictor variables, 𝜖𝜖 is the residuals derived
from the differences between tabulated values and predicted values.

Assumptions

For a linear regression model with ordinary least squared method, there are several
assumptions:
1. residuals are approximately independent, homoscedastic and normally
distributed with zero mean; and2. there is linearity between the response and predictors.

Model adequacy assessment approach
Linear least squared was used to estimate the parameters in the regression. The model
will be evaluated using the coefficient of determination
and F-test to assess the model
adequacy.
could give the strength of the linear relationship between the predictor andthe response variable. F-test provides the test result of whether the linear regression
model provides a better fit to the data than a model only with intercept. The F-statistics
is:

,

where

the

, SSE(F) is the predicted full model, and SEE(R)is the reduced model with only the intercept term, df(Residual) = n – (k+1) where k isnumber of parameters being estimates and n is the sample size. The
where

the

, is the residual, is the total sum of squares.

sum

of

squares,

,
and

the

Four different types of diagnostic plots will assess each assumption of a linear regression
model:BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

102 | NEA/CSNI/R(2021)10
1. Residuals vs. fitted: is used to check the residual homoscedasticity assumption.
If a linear model is correctly specified, the residual vs. fitted plots should nothave any systematic features with a flat line. Studentised Breusch-Pagan tests
would be used if this graph does not show a visible flat line with separated dots
around it. The is the studentised test statistics and it is calculated by findingout thewhen running the auxiliary regression equation of; examined by regressing the squared residuals on eachindependent variable.
2. Normal QQ: is used to check whether the residuals are normally distributed.
A Shapiro-Wilk normality test is another option used to check this assumption.
The W statistic is:

Where theare the ordered sample values, are constants generated fromthe means, variances and covariance’s of the order statistics of the sample froma normal distribution.
3. Residuals vs. leverage: is used to check the outliers and identify the influential
cases, since extreme values may affect the regression results if those values are
excluded or included from the analysis. Outliers could be identified if they areoutside of the Cook’s distance dashed line.

Case 1
The hazard frequency/magnitude model here for Case 1 fits the linear regression modelwith log transformation of return period to describe the relationship between magnitude
and return period since the synthetic model provided is linear. The difference is theresidual term in regression equation, which is the vector values of the differences
between observed values and predicted values� 𝚤𝚤 = 𝛽𝛽�0 + 𝛽𝛽�1 𝑙𝑙𝑙𝑙𝑙𝑙(𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝1 ) + 𝜖𝜖𝑖𝑖
𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀

Model adequacy assessment approach
The estimated equation is:.
The estimated values above are the least square estimates of the intercept and slope. They
have standard error for the intercept and slope are 0.07006 and 0.03083, respectively.
The coefficient of determination, which means 99.88% of theuncertainty in the predicted values can be explained by the straight-line regression
model. The assessment of model uncertainties is essential to the practical application of
the models to inform the decision makers. The small uncertainties of the models will
provide the confidence on decision making.
Figure D.1 shows the tabulated/predicted data and the regression model. There is anarrow range of 95% confidence band. Two tabulated values out of ten not underlying
on the 95% confidence band, but they all fall into the prediction intervals.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 103

Figure D.1. Return period versus mean magnitude

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

104 | NEA/CSNI/R(2021)10

The following, Figure D.2 to Figure D.5, are the four different diagnostic plots for the
regression model:
Figure D.2 is the plot of residuals vs. fitted values. The residuals do not change a lot as
the fitted values increase or decrease. Here, an approximate horizontal line without any
patterns shown in Figure D.2 indicates a linear relationship. Figure D.3 is the normal QQ
plot. Although there are a few points that do not follow the straight dashed line, thenormality of assumption is not violated. A Shapiro-Wilk normality test would be used
to check this assumption. p-value (0.161) of this test suggests that no evidence rejected
the normality assumption. Figure D.4 is the plot of Scale-Location.
An approximate flat line shown in this figure demonstrates a constant variance of theresiduals. Figure D.5 is the residual vs. leverage plot. There is a data point in this figurethat exceeds Cook’s distance, which is probably influential to the regression results.
Excluding this outlier data point and re-running the log-transformed variable model, we
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 105

see that there is not a huge difference between the two models. The sample size is small
and there is a large time interval corresponding to each observation. The return period
for expected magnitude prediction on this benchmark is larger than 10 000, which means
excluding this point is not appropriate.
Figure D.2. Residuals vs fitted

Figure D.3. Normal Q-Q

Figure D.4. Scale-location

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

106 | NEA/CSNI/R(2021)10
Figure D.5. Residuals vs leverage

Results
The following table shows the magnitude prediction for Case 1 for the return period of:Table D.1. Resulting magnitude prediction for Case 1
Return Period (years)
Magnitude (metres) Case 1

500
1.87

5 000
2.38

50 0002.88

500 000
3.39

The model for Case 1 indicates that for 1% increase in the return period variable, themagnitude increases by about 0.506*log (1.01) = 0.005.

Results of Assessment
Utilising the linear regression to predict magnitude is simple and efficient, there are sometechnical basis and limitations. For example, the given range of magnitudes tabulated for
up to 10 000 return period years and limited observations cause high confidence andmodel uncertainty. If a linear regression is correctly specified, all of its assumptions must
be met including the linearity of data, normality of residuals, independence andhomogeneity of residuals. For most of the cases, the data does not simply follow thelinearity assumption, we need to transform their forms to meet the assumptions. Besides,
even though our model for Case 1 has high value of , the small size and outliers
happened may influence the model results so that there are a few points not underlying
on the confidence intervals. The biggest difference between the synthetic model and theregression we fit is the random error term associated with the model, which makes theregression model be a stochastic model. However, the trend of the linear regression isconstant, so that an exact relationship between variables is determined without
considering the random error component. The linear regression model hypothesises a
probabilistic relationship between magnitude and return period.

Case 2a
The hazard frequency/magnitude model here for Case 2 does not have a synthetic model
provided. For the first case, Non-linear regression would be considered to fit themagnitude vs. return period relationship for Case 2a. A non-linear regression model has
the following general form:
�x) + 𝜖𝜖𝑖𝑖
𝑌𝑌𝑖𝑖 = 𝑓𝑓(θ,BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 107

Where Y is the response variable, x is the independent variable, θ is the vector of model
parameterthe relationship between x and y through the function 𝑓𝑓(θ, x). and 𝜖𝜖 is theresidual error term.

Case 2a has an exponential form relationship between outcome interest and predictor,
the estimated equation is:� 𝚤𝚤 = −36.9𝑒𝑒 −0.000009781∗𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑖𝑖 + 37.46
𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀

Figure D.6 shows the tabulated/predicted data and the regression model with 95%
confidence interval and prediction interval. There is a narrow range of 95% confidenceband when the return period is below 1 000 years, and there are a few tabulated values
do not fall on the 95% confidence band.
Figure D.6. Magnitude vs return period

Assumptions
For a non-linear regression model with least squared method also has the same residuals
assumption:
Residuals are approximately independent, homoscedastic and normally distributed withzero mean.

Model adequacy assessment approach
Unlike the linear regression, the 𝑅𝑅 2 could explain the goodness of fit for a model, 𝑅𝑅 2 isinvalid for non-linear regression. To access the goodness of fit for non-linear regression,
this could be done by looking that the correlation between the actual observed values

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

108 | NEA/CSNI/R(2021)10
and values predicted by a model. A high correlation coefficient (0.99) indicates thatobserved and predicted values are very close to each other.
Since the non-linear regression has the same residuals assumption as the linear
regression has, the following Residuals vs. Fitted, normal Q-Q plot, Autocorrelation
plots would check the residual assumption separately.
Figure D.7. Residuals

Figure D.8. Autocorrelation

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 109

Figure D.9. Normal Q-Q plot of standardised residuals

For the Residual vs Fitted plot, each corresponding residuals from fitted values does notchange a lot between each other, which indicates the assumption of residualhomoscedastic is met. However, it seems like the autocorrelation plot show an increasing
trend. A Runs test was used to determine the autocorrelation exists. The test statistics of𝑟𝑟− 𝜇𝜇𝑟𝑟
2𝑛𝑛 𝑛𝑛Runs test is 𝑍𝑍 =, where 𝜇𝜇𝑟𝑟 = 𝑁𝑁0 1 + 1 is the expected number of runs, 𝜎𝜎𝑟𝑟 =𝜎𝜎𝑟𝑟

2𝑛𝑛0 𝑛𝑛1 (2𝑛𝑛0 𝑛𝑛1 −𝑁𝑁)is the standard deviation of the number of runs, r is the observed number
𝑁𝑁 2 (𝑁𝑁−1)

of runs, 𝑛𝑛0 is the number of observations below the threshold, 𝑛𝑛1 is the number ofobservations above the threshold, and 𝑁𝑁 = 𝑛𝑛0 + 𝑛𝑛1 . A p-value of 0.18 indicatesautocorrelation assumption is not violated. The last graph is normal Q-Q plot with a pvalue (0.27) of Shapiro-Wilk normality test indicates a normality assumption is met.

Results
The following table shows the magnitude prediction for Case 2-1 for the return periodof:Table D.2. Resulting magnitude prediction for Case 2-1
Return Period (years)
Magnitude (metres) Case 1

500
0.74

5 000
2.32

50 00014.83

500 000
37.18

Results of assessment
An exponential non-linear regression is a good option to fit the non-linear relationshipbetween predictors and outcome and to meet the desired condition that the magnitude
would be stabilised when the return period is large enough. However, this technique may
have some technique limitations for Case 2a. First of all, the limited observed values andtheir trend corresponding to the return period is dramatically increasing, which makes it
hard to estimate the maximum mean magnitude when the return period is large enough.
Due to the small sample size, a non-significant p-value of runs test could not detect theresiduals autocorrelation issues for our non-linear regression model. The power for our
residual correlation is only 0.52. If we want to achieve our power to 0.8, the sample sizeof residuals must be at least 22. Besides, the confidence and prediction intervals for a
return period after 2 500 years has wider bands than those intervals for return period
before 2 500 years. The small sample size with a range of magnitudes is tabulated fromBENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

110 | NEA/CSNI/R(2021)10
long-term return periods, which demonstrates a larger uncertainty of a model. In
addition, compared with the log-transformed linear regression for Case 1, exponential
non-linear regression has larger predictions for the long-term return period, but has
similar predictions for the short-term return period.

Case 2b
For Case 2b, a simple linear regression with the Cochrane-Orcutt Procedure would beconsidered. The estimated equation is:
𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑖𝑖 = 0.586522 + 0.000339738(𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟𝑟 𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑖𝑖 )The Cochrane-Orcutt Procedure could adjust the serial correlation of errors for a linear
regression model. The residuals generated from the simple linear regression is first-order
autoregressive structure based on the following two graphs, which could be written as𝜖𝜖𝑡𝑡 = 𝜌𝜌𝜖𝜖𝑡𝑡−1 + 𝜖𝜖𝑡𝑡 , |𝜌𝜌| < 1 , and the model could be transformed as 𝑦𝑦𝑡𝑡 − 𝜌𝜌𝑦𝑦𝑡𝑡−1 =𝛼𝛼(1 − 𝜌𝜌) + 𝛽𝛽(𝑋𝑋𝑡𝑡 − 𝜌𝜌𝑋𝑋𝑡𝑡−1 ) + 𝜖𝜖𝑡𝑡 .
Figure D.10. Series redis(mod) [1:26]

Figure D.11. Series redis(mod) [1:26]

The coefficient of determination is 0.9994, which means 99.94% of the uncertainty inthe predicted values can be explained by the linear regression model. Figure D.12 shows
the tabulated/predicted data and the regression model. It seems like a few points do not
fall into the 95% confidence band.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 111

Figure D.12. Return period vs mean magnitude

The following figures are diagnostic plots for the regression model.
The first plot is residuals vs. fitted values. In this plot, it seems like the line has somepatterns shown below. However, each of corresponding residuals from each fitted value
only changes a very small range, which means no large variance between the fitted
values and tabulated values. The second graph, the normal Q-Q plot, noticed that thefirst and the last residuals do not follow the referenced dashed line. The last figure ischecking the residuals correlation assumption. There is no trend shown on this graph.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

112 | NEA/CSNI/R(2021)10
Figure D.13. Residuals versus fitted

Figure D.14. Normal Q-Q

Figure D.15. Autocorrelation

Results
The following table shows the magnitude prediction for Case 2b for the return period of:Table D.3. Resulting magnitude prediction for Case 2b
Return Period (years)
Magnitude (metres)

500
0.75

5 000
2.28

50 00017.57

500 000
170.46BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 113

Results of assessment
Although this simple regression with the Cochrane-Orcutt Procedure solved the residuals
correlation issue, the prediction of magnitude has larger differences for long-term
periods compared with other cases, and the normality assumption is violated due to thefirst and the last residuals points. If the sample size is small, the Type I error rates will
not be far from the target significance level, but the non-normality residuals would
generate the inaccurate prediction interval. Besides, the model could not obtain an ideal
and stabilised estimated magnitude for a long time period. If we could have more
observations, especially for long time-period, we could probably estimate a model for
ideal conditions.

Case 2c
The hazard frequency/magnitude model here for Case 2c has the estimation of theuncertainty on the magnitude for long time intervals. However, the uncertainty of shorttime intervals has not been provided for Case 2c. Therefore, the first step is using simple
linear regression to simulate those unknown uncertainties for short time intervals. Basedon the known mean magnitude, the standard deviation, 5th percentile and 95th percentile
values for long time intervals used three different linear regression models to estimate
the coefficient of relationship between the mean magnitude and those uncertainty
characteristics. A complete table is shown below:
Table D.4. Benchmark data provided for Case 2-3
Return Period (years)
Magnitude (metres)

Mean

1
053

2
0.53

5
0.54

10
0.55

50
0.59

100
0.62

500
0.79

1 000
0.95

3 000
1.6

10 0004.0

SDev

0.00670.006770.008080.009390.014620.01854

0.04

0.06

0.15

0.46

5th

0.508

0.50820.51600.52370.5546

0.57787

0.72

0.85

1.3

3.2

95th

0.5564

0.55640.56830.58030.62790.66376

0.85

1.1

1.8

4.7

To create the predictive model based on the return period and the uncertainty of meanmagnitude characteristics, the standard deviation, and the 5th and 95th percentile values,
the multicollinearity of those predictors is a big issue. A principal component analysis is
to find a new linear combination of those variables that contains much information bylooking at the highest variance. The new variables derived by the principal component
analysis are orthogonal, where the first axis contains the most information, and thesecond axis contains the second most information, and so forth. The following plotshows that only one component explains around 93% variance, and two components
explain around 99% variance in the data set.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

114 | NEA/CSNI/R(2021)10
Figure D.16. Principal component analysis results

The next step is to regress the response variable mean magnitude in a set of principal
components obtained from the return period, standard deviation, and 5th and 95th
percentiles values of the mean magnitude. The principal component regression with afew components often could explain most of the variability of all predictors andconsequently with a relationship with outcome interest. The first component from theprincipal component analysis would be considered to be used in the regression model,
and a polynomial term of the first component would be added in order to meet thenormality assumption. The principal components regression model for Case 2c has thefollowing form:2�
log (𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑀𝑖𝑖 ) = 1.04862 + 1.56278𝑃𝑃𝑃𝑃1 − 0.18800𝑃𝑃𝑃𝑃1, which means 98.62% of
The coefficient of determination 𝑅𝑅 2 = 0.9862
the uncertainty in the predicted values can be explained by the squared regression model.
Besides, F-statistics (428.7) and significant p-value (<0.05) points out the model
provides better fit than the null model.
Figure D.17 shows the tabulated/predicted data and the regression model. Predictedvalues for long time periods could not fall into the 95% confidence interval but fall into
the prediction interval.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 115

Figure D.17. Return period versus mean magnitude

Since principal component regression still uses the linear regression to predict the
response variable by using the components as predictors, residuals vs. fitted and normal
Q-Q plots could be utilised to check the residuals’ homoscedasticity and normality
assumption. Another Durbin Watson test would check the measurement of
autocorrelation in residuals from regression model. The test statistic is 𝐷𝐷𝐷𝐷 =2∑𝑇𝑇
𝑡𝑡=2(𝑒𝑒𝑡𝑡 − 𝑒𝑒𝑡𝑡−1 ), where the 𝑒𝑒𝑡𝑡 are residuals. Diagnostic plots for the regression model are2∑𝑇𝑇𝑒𝑒
𝑡𝑡=1 𝑡𝑡

shown below:

Figure D.18. Residuals versus fitted

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

116 | NEA/CSNI/R(2021)10
Figure D.19. Normal Q-Q

For the residuals vs. fitted values plot, it seems like the line has some patterns shown
between fitted values and residuals. However, residual values do not change a lot
between each other, which means there are small variances between the fitted values andtabulated values. For the normal Q-Q plot, a few residuals points do not completely
follow the referenced dashed line. Using the Shapiro-Wilk test again to ensure thenormality assumption test, a p-value with 0.18 (>0.05) suggests that the null hypothesisof normality assumption could not be rejected. Besides, the Durbin Watson test with pvalue 0.09 also indicates the autocorrelation assumption is not violated.

Results
The following table shows the magnitude prediction for Case 2-3 for return periods of:Table D.5. Resulting magnitude prediction for Case 2-3
Return Period (years)
Magnitude (metres) Case 1

500
0.71

5 000
2.01

50 00013.75

500 000
37.09

Results of assessment
This principal component regression has a few technical limitations and uncertainties.
First, the uncertainty of mean magnitude for a short time interval has not been provided
in the table, even though the unknown standard deviation of each mean magnitude issmall. In order to estimate the standard deviation, the 5th and 9th percentile of meanmagnitude, an imputation has to utilised, which produced some uncertainty for those
imputed values. To co-operate that uncertainty and return period, we used the principal
component regression that only one component of those combinations represents 93%
of variances in the dataset, inducing larger uncertainty in our model. Compared with theresults estimation obtained from Case 2a, we noticed that the confidence and prediction
intervals for Case 2c are all wider than those for Case 2a. Model fitting also does not
perform well when there is a long time return period so that the predicted values do not
fall into the 95% confidence interval. Besides, this principal component regression
model could not meet the expected condition that the mean magnitude would keepincreasing as the return period increases.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 117

Annex E. Submission by IRSN
E.1.1. Case 1 – Known model producing the synthetic data
Table 1.
Synthetic data
for Case 1

t
1
2
5
10
50
100
500
1000
2000
10000

M
0.50
0.65
0.85
1.00
1.40
1.50
1.90
2.00
2.20
2.50

The synthetic model (SM-Case1):
M = 0.5 + 0.5 × log10 ( a × t )The synthetic model serves as surrogate for acomplex phenomenological process.
Inputs:Return time intervals (or periods) t ;Outputs:
Magnitude

M for an annual maxima (AM)event. Variable a equal to 1 in cas1:
M = 0.5 + 0.5 × log10 ( t ) ).

Figure 1. Hazard curve (SM-Case1).

Benchmark completion: Provide, using the data for Case 1, a model that best describes
the frequency/magnitude relationship and the associated analysis and insights.
The randomness, homogeneity and stationarity of data are necessary conditions to
conduct a frequency analysis. As a magnitude M produced by the synthetic model issupposed to be for annual maxima (AM) events (mean rate of events λ =1 , this parameter
gives the number of events per year), the hypothetical data set should be independent.
They should also be homogeneous since the magnitudes are produced by one synthetic
model (the same statistical population). The sample size is a prerequisite for a frequency
analysis, as well. This last-mentioned condition can be very easily satisfied since weknow the synthetic model producing the data (we can produce as much data asnecessary).
Assumption 1: The magnitudes are assumed to be stationary,
independent and homogeneous.
As mentioned earlier, a magnitude M is produced by the synthetic model with a rateof events equal to one ( λ =1 ). Return periods can then be estimated with one of the two
following frequency models:
1. An annual maxima (AM) frequency model in which the distribution of the AM
events converges to a GEV one.
2. A Peaks-Over-Threshold (POT) model in which the distribution of theexceedances over the threshold u converges to an exponential one (ofparameter θ ).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

118 | NEA/CSNI/R(2021)10
Let xT be the event with the return period T and F the non-exceedance probability
function. In its general definition, xT corresponds to the quantile of probability ofexceedance equal to 1 λT (i.e. equal to 1 T when λ =1 ): Pr ( X ≥ xT ) =1 T ( F isF ( xT ) =λthe survival function).
In the 1st case (SM-Case1), a set of 100 return periods is randomly sampled. Associatedmagnitudes are calculated with the synthetic model and a classic frequency estimationis then performed with a GEV distribution. To take this issue one step further, a simple
analogy of the synthetic model with GEV and GP distribution functions is performed.

1. AM/GEV frequency model:
The GEV distribution is the limiting distribution for the AM of independent andidentically distributed (iid) random variables.

 −1+ ξ x − µ − 1ξ  σ ξ ≠0
F ( x ) = e
 − e− ( x−µ ) σ
ξ =0
e

σ >0

Where µ ,and ξ are the location, scale and shape parameters, respectively. A
sample of 100 magnitudes (AM events with an effective duration equal to 100 years) is
generated and analysed. The hypothetical observational data set (from a synthetic model)
is selected in such a way that their empirical distribution looks as natural as possible(close to the natural behaviour and variability of external hazards such as floods, extreme
temperatures or high winds, etc.). More concretely, the sampling is done as follows: theR “sample” function is used to randomly sample 100 values of “return times” t from 1
to 100 years. In all, 95 values have empirical return periods less than 100 years. It is
worth noting that the case with aleatory sampling in a period of 10 000 years is notrealistic and cannot represent natural hazards.
library("evd")
library("ggplot2")# Hypothetical data from synthetic model# A sample of 100 AM in which the 100-year return level was not exceededw <- 100
tt1 <- sample(1:100, size=w, replace=TRUE);
xt <- round(0.5+0.5*(log10(tt1)),2)# Plotting positions - empirical probabilities (Weibul)
Px <- 1:length(xt)/(length(xt)+1)# associated empirical return periods
Tx <- 1/(1-Px)# Theoretical fitting
fit.GEV <- fevd(xt, type="GEV")# Estimated parameters: # location scaleshape#
1.2589728 0.1924330 -0.2943914
ttic <- c(1.25, 1.5, 2, 3, 5, 10, 20, 50, 100, 500, 1000, 2000, 3000)
RLs <- ci(fit.GEV, type = "return.level", return.period=ttic, method ="normal", alpha=0.30)

Plot the fitting . . .

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 119
data1 <- data.frame(xx = Tx, yy = sort(xt))
data2 <- data.frame(xxx = ttic, yyy = RLs [,2])data3 <- data.frame(xxxx = ttic, lb = RLs [,1])data4 <- data.frame(xxxx = ttic, ub = RLs [,3])sp <- ggplot(data.frame("pp"=NULL, "tf"=NULL))
sp + geom_point(aes(x=data1$xx, y=data1$yy, colour="PPcol" )) +geom_line(aes(data2$xxx, data2$yyy,colour="Fitcol",linetype="Fitline1")) +geom_line(aes(data3$xxxx, data3$lb, colour="ICcol", linetype="Fitline2"))+
geom_line(aes(data4$xxxx,data4$ub,colour = "ICcol",linetype="Fitline2")) +scale_color_manual("",labels=c("PP","Fit","IC70%"),breaks=c("PPcol","Fitcol","ICcol"),values=c("black","black","black")) +scale_linetype_manual("",breaks = c("fitline1","fitline2"),
values = c("solid","longdash")) +guides(fill = guide_legend(order = 1), colour = guide_legend(order = 2, keywidth =2.5, override.aes = list(linetype =c("blank", "solid", "longdash"),shape=c(16,NA,NA)))) +scale_x_log10(breaks=c(10, 100, 500, 1000), labels=as.character(c(10, 100, 500,
1000)), name = "Return period (years)") +scale_y_continuous(name = "Magnitude (m)", limits = c(0,2)) +ggtitle("Fitting with a GEV distribution (AM data)") +theme(
plot.title = element_text(color = "Black", size=10, face = "bold.italic"))

The fitting is bounded to a final value equal to µ − σ ξ (equal to 1.91 m). The 500 000year return level is equal to this end value (1.90 m) which is much lower than the 3.35 m
calculated with the synthetic model for the same return period. Paradoxically, it can beseen from Figure E.1 that the adequacy of the GEV distribution is quite good and theuncertainty (in term of confidence intervals) is very low for high return levels.
Figure E.1. Fitting with a GEV distribution with the SM-Case1 data sets (µ = 1.26; σ = 0.19; ξ = −0.29 )The GEV distribution combines three distributions, identified by Fisher and Tippet(1928) into a single form. Indeed, depending on the value of the shape parameter ξ , theGEV can take the form of the Gumbel ( ξ = 0 ), Fréchet ( ξ > 0 ) or the Reverse Weibull

distributions ( ξ < 0 ). On the other hand, and as mentioned earlier in this section, the scale

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

120 | NEA/CSNI/R(2021)10

σ > 0 ). The probability andparameter of the GEV distribution should be positive (quantile functions for the GEV distribution can be written as:
A linear relationship between magnitudes M and thelogarithm of return periods t could be obtained if the xaxis was plotted in a logarithmic scale. Note that thislinear relationship cannot be obtained with all the
asymptotic extreme value distributions and only the
exponential form (the Gumbel distribution) might besuitable (but with the limitation of bad fitting of small
magnitudes associated to small return periods up to 50100 years).

1 −1+ ξ x −µ − ξ
 F ( x ) = e  σ ξ≠0
− x −µ σ−e ( )

ξ =0

e

−ξ
σ

 xF = µ − 1 − {− log (1 − F )} ξAs can be seen from the probability equation, the GEV distribution can as well have theform of the synthetic model when the shape parameter is equal to -1. With such a value,
the theoretical upper tail can only be finite and bounded (as may be useful for estimatesof specific cases of extreme values which may have an upper bound, as is the case here).
This hypothesis is in line with the bounded theoretical fitting presented in Figure E.1.
1
xF = µ + σ 1 + log (1 − F )  → xT = µ + σ + σ log   → xT = µ + σ − σ log (T )T xT = µ + σ − σ log (10 ) × log10 (T )This last equation is similar to the proposed synthetic model as follows:

It can be concluded from these two conditions that, with a shape parameter equal to -1,
the scale parameter can only be negative. As just noted, a GEV distribution obviously
cannot be used with a negative scale parameter: then it cannot describe thefrequency/magnitude relationship for the SM-Case1.

2. The POT GPD/Exponential frequency model:
The GP distribution for the return levels (with an exponential distribution for theexceedances over the threshold u ). The exceedances over the threshold u follow an
exponential distribution of parameter ρ .

Fu ( X ) = 1 − exp  −ρ ( x − u ) 
1
1
1
 1 
xT − u =− log (1 − F ) → xT =u − log  → xT =u + log ( λT )ρρ  λT ρ

Therefore, theT − year return level can be written as:

11
xT = u + ln ( λ ) + log (10 ) log10 (T )ρρ

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 121

This last equation is similar to the proposed synthetic model as follows:

Considering the case in which the rate of events is equal to one ( λ =1 ), the first condition
gives a threshold u = 0.5 . It can also be easily concluded from the second condition that
=
ρ log
=
(10 ) 0.5 4.605 .
Assumption 2: The number of events is randomly sampled
following a Poisson Process. It is also assumed that the rate ofevents is equal to 1.
In the first step, the developed frequency model is used to describe thefrequency/magnitude relationship. Magnitudes are then sampled on a period of time w, say 100 years.
set.seed(1248)# effective durationw <- 100# Thresholdu <- 0.5# annual rate of events
lambda <- 1
# Parameter of the exponential distribution
rho <- log(10)/0.5# Assuming the Nbre of arrivals happen randomly following a Poisson Process
N <- rpois(n = 1, lambda = lambda * w)
# Return levels (exceedances follow an exponential distribution)
xt <- u + rexp(n = N, rate = rho)

In the next step, the fitting with the theoretical distribution (exponential) is performed.
The fitting of the SM-Case1 data set (w=100 years) with a GPD distribution (=u 0.5 =& λ 1 ) for the return levels (the exceedances over u are exponential) ispresented in Figure E.2. The results are also shown in Table E.1. A good and visuallyadequate fitting is obtained. Indeed, all the observed probabilities are in the 70%
confidence interval. As shown in Table E.1, the 500 000-year return level is equal to3.19 m which is close to the 3.35 m calculated by the synthetic model for the same return
period. On the other hand, the uncertainty (in term of confidence intervals) is reasonably
low for high return levels. All the simulations were carried out within the R environment(open-source software for statistical computing: http://www.r-project.org/). The Renext
library (IRSN and Alpstat, 2013), developed by the French Institute for RadiologicalProtection and Nuclear Safety, was used for frequency estimations. The Renext package
was specifically developed for flood frequency analyses using the POT method.
library("Renext")
library("ggplot2")
fit.exp <- Renouv(xt, effDuration = w,
"exponential",
Tlim= c(1, 10000),
plot = TRUE)

distname.y =threshold = u,pct.conf = 70,

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

122 | NEA/CSNI/R(2021)10
Figure E.2. Fitting of the SM-Case1 data sets (w=100 years) with GPD distribution (𝒖𝒖 = 0.5 & 𝛌𝛌 = 𝟏𝟏) for
the return levels (the exceedances over u are exponential).

Table E.1. Comparison of magnitudes (m) calculated with the synthetic data for Case 1. The values inbrackets correspond to the 70% confidence intervals.
Return Period (years)
Synthetic data for Case 1
POT model (GPD/Exp)

500
1.85
1.78 (1.66-1.92)

5 000
2.35
2.25 (2.08-2.45)

50 0002.85
2.72 (2.51-2.98)

500 000
3.35
3.19 (2.94-3.50)

3. Case 1 – Conclusion
In this first case, both the hypothetical data and the synthetic model are known. It can beconcluded that the GEV distribution is not very suitable to describe all thefrequency/magnitude relationships for the SM-Case1 data. These data are best describedwith a POT frequency model in which the distribution of the return levels converges toa GPD (threshold equal to 0.5 and rate of events equal to 1) with an exponential
behaviour of the exceedances over the threshold. Indeed, the desired high return levels
are estimated with this frequency model and compared to the hypothetical observational
data from the synthetic model. The relative difference in magnitudes did not exceed 5%(7 cm, 10 cm, 13 cm and 16 cm for the 500-, 5 000-, 50 000- and 500 000 year returnlevels, respectively). Moreover, all the plotting positions are inside the confidence
interval despite the fact that the latter is very narrow (low uncertainty).

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 123

E1.2. Case 2 – Unknown model producing the synthetic data
For Case 2, there is no synthetic model provided. Furthermore, three parts to thisexample are provided.
Case 2(a): provides the synthetic data (10 data points) with no uncertainty provided;Case 2(b): provides additional synthetic data (26 data points) with no uncertainty
provided;Case 2(c): provides the synthetic data (10 data points) with uncertainty estimates onsome of the data.
We need to provide, using the data for Case 2, a model that best describes thefrequency/magnitude relationship and the associated analysis and insights.
Non-linear least-squares estimates of the distribution parameters are performed using the“nls” R function. Confidence intervals are then calculated and plotted. The same
assumptions are used for this second case.
Assumption 1: The magnitudes are assumed to be stationary,
independent and homogeneous.
Only parts 2a and 2b are evaluated herein. As in Case 1, return periods can then be
estimated with one of the two following frequency models:
1. The AM/GEV frequency model: Annual maxima (AM) frequency model inwhich the distribution of the AM events converges to a GEV one.
2. POT/GPD frequency model: Peaks-Over-Threshold (POT) model in which the
distribution of the exceedances converge to a GPD.
For both Cases 2a and 2b, non-linear least-squares estimates of the GEV and GPD
parameters are performed. It is worth noting that the GPD must give almost the same
parameters and fitting (results not presented hereafter). The fitting with confidence
intervals is presented in Figure E.5 and Figure E.8. The adequacy of the theoretical
distribution in both Cases 2a and 2b is visually quite good with heavy tails (very highshape parameter ξ ≈ 0.96 ). A comparison of T-year return levels (corresponding to 500-,
5 000-, 50 000- and 500 000-year return periods) for both Case 2a and Case 2b arepresented in Table E.2 (the values in brackets correspond to the absolute and relative
widths of the 70% confidence intervals). The results in Table E.2 indicate that therelative widths of confidence intervals in Case 2b (with 26 synthetic data points) are
1.5 times narrower than those obtained in Case 2a (with 10 synthetic data points) for the
GEV distribution.
Table E.2. The T-year quantiles with absolute and relative widths of their 70% confidence intervals
Return(years)Period500

5 000

50 000500 000

GEV for Case 2-a0.75 (0.68-0.84)(21.3%)

2.32 (1.72-3.21)(64.2%)

16.79 (10.75-26.44)(93.4%)

150.34 (89.35-254.16)(109.6%)

GEV for Case 2-b0.76 (0.71-0.82)(14.5%)

2.31 (1.88-2.87)(42.8%)

16.50 (12.17-22.46)(62.4%)

146.53 (102.70-209.56)(73.0%)BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

124 | NEA/CSNI/R(2021)10

Cases 2a and 2b – Conclusion
In these two cases, only the synthetic data are known (no synthetic model provided). It
can be concluded that both the MA/GEV and POT/GP distributions best describe themagnitude/frequency relationship with heavy tails (very high shape parameter, thenthere is no physical limit!). On the other hand, as more data are provided in Case 2b, theconfidence intervals are narrower (reduced uncertainty).
It is obvious that, with such a high shape parameter, these data reflect what can beobserved for natural hazards up to a return period of 1 000 - 10 000 years. Beyond this,return magnitudes increase much more quickly to very high values. This is the main
characteristic of models with very heavy tails (as is the case here).

Case 2a: the MA/GEV frequency model:
10 synthetic data points are provided with no uncertainty provided;# Benchmark on External Events Hazard Frequency and Magnitude #Statistical Modelling# Case 2 – Unknown Model Producing the Synthetic Datalibrary("evd")
library("ggplot2")# Case 2-a synthetic datat <-c(1,2,5,10,50,100,500,1000,3000,10000)
M <- c(0.53,0.53,0.54,0.55,0.59,0.62,0.79,0.95,
1.60,4.00)# Plot data - Hazards curve
qplot(t, M, geom=c("point", "line"), xlim = c(0, 10000), ylim = c(0.5, 4),main = "Data (with no synthetic model provided) - Cas2-a", xlab ="Return period (years)",ylab = "Magnitude (m)")

Table E.3. Synthetic data for Case 2a
t
1
2
5
10
50
100
500
1 000
3 000
10 000

M
0.53
0.53
0.54
0.55
0.59
0.62
0.79
0.95
1.57
3.97

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 125

Figure E.3. Hazard curve (Case 2a synthetic data)

Figure E.4. Fitting Case 2a synthetic data

df1 <- data.frame(t = t, xt = M)
# Fitting with a Non-linear least-squares estimates of the GEV parameters
fit <- nls(xt ~ qgev(1-1/t, loc = mu, scale = sigma, shape = xi), data=df1[-1,], start= list(mu = 0.12, sigma = 0.2, xi = 1))
coef(fit)
# mu
sigma
xi
# 0.5593335970 0.0004576594 0.9649876346# Plot fitting (code not presented herein)

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

126 | NEA/CSNI/R(2021)10==
σ 0.0005;=
ξ 0.9650 ))Figure E.5. Fitting Case 2a synthetic data (GEV
( µ 0.5593;

Case 2b: the MA/GEV frequency model:
26 synthetic data points are provided with no uncertainty provided;
library("evd")
library("ggplot2")# Case 2-b synthetic datat <-c(1,2,5,10,15,20, . . . ,750,1000,3000,10000)
M <- c(0.53,0.53,0.54,0.55, . . . ,0.95,1.57,3.97)# Plot data - Hazards curve
qplot(t, M, geom=c("point", "line"), xlim = c(0, 10000), ylim =c(0.5, 4), main = "Data (with no synthetic model provided)- Cas2-b", xlab = "Return period (years)",ylab ="Magnitude (m)")

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 127

Table E.4. Synthetic data for Case 2b
t
1
2
5
10
15
20
25

M
0.53
0.53
0.54
0.55
0.56
0.56
0.57

30
40
50
60
70
80
90
100
125
150
175
200
300
400
500
750
1 000
3 000
10 000

0.57
0.58
0.59
0.60
0.60
0.61
0.62
0.62
0.63
0.65
0.66
0.67
0.71
0.75
0.79
0.87
0.95
1.57
3.97

Figure E.6. Hazard curve (Case 2b synthetic data)

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

128 | NEA/CSNI/R(2021)10
Figure E.7. Fitting the Case 2b synthetic data

df1 <- data.frame(t = t, xt = M)
# Fittingfit <- nls(xt ~ qgev(1 - 1/t, loc = mu, scale = sigma, shape = xi),
data = df1[-1, ],start = list(mu = 0.12, sigma = 0.2, xi = 1))
coef(fit)
# mu
sigma
xi
# 0.5707180614 0.0004624945 0.9619867644# Plot fitting (code ggplot not presented herein)==
σ 0.0005;=
ξ 0.9619 )Figure E.8. Fitting Case 2a synthetic data with a GEV distribution( µ 0.5709;

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 129

Annex F. Submission by KAERIIntroduction
This benchmark study aims to apply statistical modelling for frequency and magnitude
estimation based on data for external event hazard assessment. Based on the results of thisstudy, it is believed that an approach to the quantification of external event IEs can beformulated and evaluated through the application of an effective statistical model. In thisstudy, analysis was based on two cases that considered benchmarks provided by the OECD
NEA. Each case was given a magnitude according to the return period. Based on this data,
an appropriate statistical model was applied through regression analysis for each case.
Based on the results, the magnitudes of 500, 5 000, 50 000, and 500 000 years werepredicted and presented.

Synthetic data analysis•

Study for Case 1
Table F.1 shows the synthetic data for Case 1 as provided by the OECD NEA. Using
regression analysis, log regression showed appropriate fitting results for the relationshipbetween magnitude and return period, as shown in Figure F.1.
Table F.1. Synthetic data for Case 1

Return period (years)
Original M (metres)

1
0.50

2
0.65

5
0.85

10
1.00

50
1.40

100
1.50

500
1.90

1 000
2.00

2 000
2.20

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

10 0002.50

130 | NEA/CSNI/R(2021)10
Figure F.1. Log regression fitting

As a result of the regression analysis on the magnitude of Case 1, the log regression
equation including variables A and B is shown in Eq. (1). The magnitude of the return
period from 1 to 10 000 years was estimated by Eq. (1); Table E.2 compares the values tothe original magnitude values proposed by the OECD NEA.
Case 1: 𝑀𝑀 = 0.2199 ∗ 𝑙𝑙𝑙𝑙 (𝑥𝑥) + 0.5034 (𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵) (1)Table F.2. Regression results for Case 1

Return period (years)
Original M (metres)
Log Magnitude (metres)

1
0.50
0.503

2
0.65
0.656

5
0.85
0.857

10
1.00
1.010

50
1.40
1.364

100
1.50
1.516

500
1.90
1.870

1 000
2.00
2.022

2 000
2.20
2.175

10 0002.50
2.529

An error analysis was then performed using the SUMXMY2 function to verify thestatistical justification of the estimated magnitudes. The SUMXMY2 function squares andsums the difference between two corresponding values, and therefore, the closer to 0, thesmaller the error between the two variables, and the more statistically valid the estimation
can be considered. The SUMXMY2 function is expressed as Eq. (2).
𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆𝑆2 = ∑(𝑥𝑥 − 𝑦𝑦)2(2)Here, 𝑥𝑥 is the value of the original magnitude, and 𝑦𝑦 is the value of the estimatedmagnitude. As a result of error analysis using the SUMXMY2 function, the sum squareerror (SSE) value was calculated as 0.005, which is very close to zero. It was therefore
judged that the estimated magnitude values were very similar to the original values andvalid. However, to minimise SSE and more precisely estimate the magnitude values, asolver function was used. The target of the SSE value was set to 0, and an optimisation
analysis was performed on parameters A and B of Eq. (1). The results are shown in Table
F.3.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 131

Table F.3. Optimisation for parameters
Parameter
A
B

Original
0.21990.5034

SSE_Solver
0.219861834
0.503363549

These optimised parameters are then used in the regression equation for Case 1, as shownin Eq. (3) below. Table F.4 compares the magnitude values estimated by Eq. (3) with thoseestimated by Eq. (1) and the original magnitude values.
Case 1: 𝑀𝑀 = 0.219861834 ∗ 𝑙𝑙𝑙𝑙 (𝑥𝑥) + 0.503363549 (𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵) (3)Table F.4. Comparison of magnitudes

Return period (years)
Original M (metres)
Log Magnitude (metres)
Optimisedmagnitude (metres)

1
0.500
0.503
0.503

2
0.650
0.656
0.656

5
0.850
0.857
0.857

10
1.000
1.010
1.010

50
1.400
1.364
1.363

100
1.500
1.516
1.516

500
1.900
1.870
1.870

1 000
2.000
2.022
2.022

2 000
2.200
2.175
2.175

For Case 1, the optimised parameters A and B (Table F.3) were similar to the existingvalues. Likewise, the SSE value of 0.0049 was also similar to the existing value of 0.005.
Fitting was then performed based on the optimised magnitude values; results are shown in
Figure F.2.
Figure F.2. Fitting result for Case 1

F.1.1. Case 1 result assessment
The results of the optimised fit for Case 1 were estimated to be similar to the size values
presented by the OECD NEA. The trend lines calculated from the estimated magnitude
values (Figure F.1 and Figure F.2) were also similarly estimated. Therefore, it can be
judged that the fitting result for Case 1 in this study is valid. Additionally, based on the
BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

10 000
2.500
2.529
2.528

132 | NEA/CSNI/R(2021)10
regression analysis estimates in Case 1, magnitude values were predicted for 500, 5 000,
50 000, and 500 000 years return period. The results are presented in Table F.5 and Figure
F.3.
Table F.5. Magnitude prediction by return period (Case 1)
Return period (years)
Magnitude (metres) Case 1 Exact
Log KAERI mean (metres)
Optimised KAERI mean (metres)

500
1.9
1.870
1.870

5 000
2.4
2.376
2.376

50 000
2.9
2.883
2.882

500 000
3.4
3.389
3.388

Figure F.3. Magnitude prediction fitting for Case 1

F.1.2. Study for Case 2
Table F.6 shows the synthetic data for Case 2 provided by the OECD NEA. In this case,regression analysis found linear regression to give the best fit between magnitude andreturn period, as shown in Figure F.4.
Table F.6. Synthetic data for Case 2
Return Period (years)
Original M (metres)

1
0.53

2
0.53

5
0.54

10
0.55

50
0.59

100
0.62

500
0.79

1 000
0.95

3 000
1.60

10 000
4.00

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 133

Figure F.4. Linear regression fitting

Equation (4) below is the linear regression equation including variables A and B as a resultof the regression analysis on the magnitude of Case 2. Magnitudes of return periods from
1 to 10 000 years were estimated by Eq. (4) and compared to the values of the original
magnitude proposed by the OECD NEA (Table F.7). In addition, 95% and 5% confidence
interval magnitude provided by the OECD NEA are also given in Table F.7, with fitting
results plotted in Figure F.5.
Case 2: 𝑀𝑀 = 0.0003 ∗ 𝑥𝑥 + 0.5651 (𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵)(4)

Table F.7. Regression results for Case 2

Return period (years)
Original M (metres)
Linear Magnitude (metres)
Original M 95% (metres)
Original M 5% (metres)

1
0.53
0.565
—
—

2
0.53
0.566
—
—

5
0.54
0.567
—
—

10
0.55
0.568
—
—

50
0.59
0.580
—
—

100
0.62
0.595
—
—

500
0.79
0.715
0.85
0.72

1 000
0.95
0.865
1.1
0.85

3 000
1.60
1.465
1.8
1.3

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

10 000
4.00
3.565
4.7
3.2

134 | NEA/CSNI/R(2021)10
Figure F.5. Linear regression fitting according to confidence interval

As a result of the regression analysis, the estimated magnitude values fell within theconfidence interval and cannot be judged as inappropriate. The estimated magnitude forthe initial return period was similar to those proposed by the OECD NEA. However, as thereturn period increases, errors in the magnitude values were found to occur. Error analysis
in this case using the SUMXMY2 function gave an SSE value of 0.22443. In order tominimise this error in Case 2, the solver function was used, where again the target of theSSE value was set to 0 and an optimisation analysis was performed on parameters A and B
from Eq. (4). The results are shown in Table F.8.
Table F.8. Optimisation for parameters
Parameter
A
B

Original
0.00030.5651

SSE_Solver
0.000344252
0.565051348

Based on the optimised parameters, linear regression analysis for Case 2 was re-estimatedvia Eq. (5). Table F.9 compares the magnitude values from the three sources [original, Eq.
(4), and Eq. (5)] along with the confidence intervals.
Case 1: 𝑀𝑀 = 0.000344252 ∗ 𝑥𝑥 + 0.565051348 (𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃𝑃 𝐴𝐴 𝑎𝑎𝑎𝑎𝑎𝑎 𝐵𝐵)(5)

Table F.9. Compare for magnitude

Return period (years)
Original M (metres)
LinearMagnitude (metres)
Optimisedmagnitude (metres)
Original M 95% (metres)
Original M 5% (metres)

1
0.53
0.565

2
0.53
0.566

5
0.54
0.567

10
0.55
0.568

50
0.59
0.58

100
0.62
0.595

500
0.79
0.715

1 000
0.95
0.865

2 000
1.6
1.465

10 000
4
3.565

0.565

0.566

0.567

0.568

0.582

0.599

0.737

0.909

1.598

4.008

—
—

—
—

—
—

—
—

—
—

—
—

0.85
0.72

1.1
0.85

1.8
1.3

4.7
3.2

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLINGNEA/CSNI/R(2021)10 | 135

In Case 2, the SSE of the optimised model was calculated to be 0.01. Accordingly, it was
judged that the parameters A and B were significantly improved compared to the existing
parameters. Fitting was then performed based on the optimised magnitude values. The
results are shown in Figure F.6.
Figure F.6. Fitting result for Case 2

F.1.3. Case 2 result assessment
It was found that when the optimisation technique was applied to Case 2, the model
performance further improved, as seen in Figure F.5. In other words, the optimised fitting
was able to estimate values similar to the magnitude values proposed by the OECD NEA.
Comparing Figure F.4 and Figure F.5, the trend lines calculated from the estimatedmagnitude values were also similarly estimated. Therefore, like Case 1, it can be judgedthat the fitting result for Case 2 in this study is valid.
Then, based on the regression analysis estimates in Case 2, magnitude values for return
periods of 500, 5 000, 50 000, and 500 000 years were predicted. The results are shown in
Table F.10 and Figure F.7.
Table F.10. Magnitude prediction by return period (Case 2)
Return period (years)
Magnitude (metres) Case 2 Exact
Linear KAERI mean (metres)
Optimised KAERI mean (metres)

500
0.78
0.715
0.737

5 000
2.2
2.065
2.286

50 00028
15.56517.778500 000
2 000
150.565
172.691

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING

136 | NEA/CSNI/R(2021)10
Figure F.7. Magnitude prediction fitting for Case 2

Conclusion
In this study, statistical analysis was applied to the estimation of two cases presented by
the OECD. In any statistical analysis, it is important to understand the characteristics of thedata set. For the given problems here, the range of the return period was 10–10 000 years,
while that of the magnitude was only 0.4–5.0 metres. Therefore, the coefficient of thesynthetic model had a great influence on the analysis results. This study demonstrates that
employing the full extent of the significant figures is important to handle the different
ranges of data values. In the future, it is expected that data-based statistical values can bebetter estimated through various verified statistical models.

BENCHMARK ON EXTERNAL EVENTS HAZARD FREQUENCY AND MAGNITUDE STATISTICAL MODELLING