\documentclass[a4paper, 12pt]{article}
\usepackage[utf8]{inputenc}
\usepackage[czech]{babel}
\usepackage{csquotes}
% Seems not necessary as we do not use qutations, but biblatex complains without it
\usepackage[a4paper, ,left=3cm,right=3cm,top=3cm,bottom=3.5cm]{geometry}
\usepackage[T1]{fontenc}      % T1 kodovani fontu pro babel cestinu
\usepackage{pdfpages} 
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{enumerate}
\usepackage{epic}
\usepackage{stmaryrd} % required for \jmp
\usepackage{subcaption}
%\usepackage{algorithm}
\usepackage{float}
\usepackage[unicode, colorlinks=true, citecolor=blue]{hyperref}
\usepackage[ruled,linesnumbered,noend]{algorithm2e}
\usepackage{bm} %bold math font
\usepackage{stmaryrd} % required for \jmp
\usepackage{tikz}
\usepackage{booktabs, multirow}
% \usepackage{acronym}
\usepackage[acronym]{glossaries}

%\usepackage{natbib}
%\bibliographystyle{abbrvnat}
%\bibliographystyle{plainnat}
%\setcitestyle{authoryear}

% TRYING TO USE BIBLATEX AND BIBER to get correct reference order
\usepackage[
    backend=biber, %natbib,
    % IMPORTANT: load a style suitable for your discipline
    %style=authoryear
    %style=authoryear-icomp
    style=bwl-FU,
]{biblatex}
\addbibresource{references.bib}% Syntax for version >= 1.2
\AtEveryBibitem{%
  \clearfield{month}%
  \clearfield{note}%
  \ifboolexpr{ 
    test {\ifentrytype{book}}
    or test {\ifentrytype{article}}
    or test {\ifentrytype{collection}}
    or test {\ifentrytype{incollection}} 
    or test {\ifentrytype{inproceedings}} 
  }{%
    \clearfield{url}%
    \clearfield{urlyear}
  }{}
}
\DeclareFieldFormat{doi}{%
  \newline
  \mkbibacro{DOI}\addcolon\space
    \ifhyperref
      {\href{http://dx.doi.org/#1}{\nolinkurl{#1}}}%
      {\nolinkurl{#1}}}
\DeclareFieldFormat{url}{%
  \newline
  \mkbibacro{URL}\addcolon\space
   \href{#1}{\nolinkurl{#1}}%   
   }

\def\Real{\mathbb R}
\def\E{{\mathsf E}}
\def\avg#1{\langle#1\rangle}
\def\abs#1{\lvert#1\rvert}
\def\prtl{\partial}
\def\unit#1{\,\mathrm{#1}}
\def\bunit#1{[\mathrm{#1}]}
\def\eps{\varepsilon}
\def\T{\intercal}
\def\grad{\nabla}
\def\div{\operatorname{div}}
\def\tr{\operatorname{tr}}
\def\vc#1{\mathbf{\boldsymbol{#1}}}     % vector
\def\tn#1{{\mathbf{\boldsymbol{#1}}}}    % tensor
\def\sqr#1{{{\left(#1\right)}^2}}
\def\todo#1{{\color{brown} TODO: #1}}
%\def\jb#1{\color{blue}#1}
\def\jb#1{{\color{green} #1}}
\def\cT{\mathcal T}
\def\cS{\mathcal S}


\DeclareMathOperator*{\argmin}{arg\,min}

\newcommand*\rot{\rotatebox{60}}

\newcommand{\js}[1]{{\color{blue}#1}} % Stebel
\newcommand{\tl}[1]{{\color{red}#1}}  % Luber
\newcommand{\jk}[1]{{\color{violet}#1}} % Kruzik
\newcommand{\pe}[1]{{\color{purple}#1}} % Exner
%\newcommand{\sy}[1]{{\color{red}#1}}  % Sysala
\newcommand{\sy}[1]{#1}  % Sysala

% plus and minus in circle
\newcommand{\opm}{
  {\mathbin{
    \mathchoice
      {\buildcirclepm{\displaystyle}{0.14ex}{0.95}{0.05ex}{.7}}
      {\buildcirclepm{\textstyle}{0.14ex}{0.95}{0.05ex}{.7}}
      {\buildcirclepm{\scriptstyle}{0.13ex}{0.955}{0.04ex}{.55}}
      {\buildcirclepm{\scriptscriptstyle}{0.08ex}{0.95}{0.03ex}{.45}}
  }} 
}
\newcommand\buildcirclepm[5]{%
  \begin{tikzpicture}[baseline=(X.base), inner sep=-#5, outer sep=-.65]
    \node[draw,circle,line width=#4] (X)  {\footnotesize\raisebox{#2}{\scalebox{#3}{$#1\pm$}}};
  \end{tikzpicture}%
}


\newcommand{\adiv}{\widetilde\div}
\newcommand{\agrad}{\widetilde\nabla}
\newcommand{\fgrad}{\overline{\grad}} 
\newcommand{\fdiv}{\overline{\div}}
\newcommand{\vnu}{\vc \nu}
\newcommand{\avgf}[1]{\{\mskip-5mu\{#1\}\mskip-5mu\}}
\newcommand{\eq}[1]{\begin{equation}#1\end{equation}}
\newcommand{\jmp}[1]{\llbracket #1 \rrbracket} % requires package stmaryrd
\newcommand{\uu}{\vc u}
%\newcommand{\alert}[1]{{\color{red}#1}}
\newcommand{\alert}[1]{#1}
\newcommand{\mbf}[1]{\mbox{\boldmath$#1$}}
\newcommand{\figpath}{graphics/}
\newcommand{\respath}{results/}
\newcommand{\figref}[1]{Obr.~\ref{#1}}
\newcommand{\degree}{^\circ}
\newcommand{\strain}{\mathbf{\boldsymbol{\epsilon}}}
\newcommand{\perm}{\mathbf{\boldsymbol{\kappa}}}
\newcommand{\permsc}{\kappa} % scalar permeability

%\hoffset=-29mm
%\textwidth=162mm
%\textheight=231mm
%\voffset=-20mm

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\title{Necertifikovaná metodika: Posouzení vlivu EDZ na transport radioaktivních látek a bezpečnost hlubinného úložiště radioaktivního odpadu pomocí výpočetních metod}
%\author{Jan Březina, Stanislav Sysala, ...}
%\date{2022}

\makeglossaries
\newcommand{\acro}[2]{\newacronym{#1}{#1}{#2}}
\acro{DAMH}{Delayed acceptance Metropolis-Hastings, algoritmus MH se zpožděným přijetím}
\acro{DFM}{Discrete fracture-matrix, kombinace kontinuálního popisu horniny a sítě diskrétních puklin}
\acro{DGM}{Discontinuous Galerkin method, nespojitá Galerkinova metoda}
\acro{EDZ}{Excavation damage zone, zóna poškození způsobená ražbou}
\acro{FETI}{Finite element tearing and interconnecting, výpočetní metoda rozložení oblasti}
\acro{GSI}{Geologický pevnostní index}
\acro{HF}{Hydrofracking}
\acro{HM}{Hydromechanika, hydromechanický}
\acro{HÚ}{Hlubinné úložiště radioaktivních odpadů}
\acro{IBUP}{Indikátor bezpečnosti úložné pozice}
\acro{MH}{Metropolis-Hastings algoritmus}
\acro{MLMC}{Multilevel Monte Carlo, víceúrovňová metoda Monte Carlo}
\acro{MPI}{Message parsing interface, knihovna pro podporu paralelních výpočtů}
\acro{MVP}{Microvelocity probe}
\acro{PVP Bukov}{Podzemní výzkumné pracoviště Bukov}
\acro{QP}{Quadratic programming, kvadratické programování}
\acro{SEPPI}{Metoda měření propustnosti hornin}
\acro{SKB}{Svensk Kärnbränslehantering, Švédská společnost pro správu radioaktivních odpadů}
\acro{SÚRAO}{Správa úložišť radioaktivních odpadů}
\acro{SW}{Software}
\acro{TBM}{Tunnel boring machine, razící štít}
\acro{THMC}{Termo-hydro-mechano-chemický proces}
\acro{TSX}{Tunnel Sealing Experiment, experiment realizovaný v Underground Research Laboratory v 
Kanadě}
\acro{URL}{Underground Research Laboratory, Kanada}
\acro{ÚOS}{Úložný obalový soubor}



\begin{document}

%\ULCornerWallPaper{1}{\includegraphics{TACR_logo.png}}
%\LLCornerWallPaper{1}{bar}
%\lipsum[1-3]

\begin{titlepage}
    \includegraphics[width=2cm]{logo_TACR_zakl.pdf}

    \vspace{2cm}
    {      
      {\huge 
      \noindent 
      oponovaná metodika \\
      \vspace{2ex}
      
      {\bf
      \noindent Posouzení vlivu\\
      zóny ovlivněné ražbou\\
      na bezpečnost hlubinného úložiště\\
      pomocí výpočetních metod}
      \par}
    }
	
    \vspace{2cm}	
    {\LARGE
      \vspace{2ex}
      {\noindent autorský tým:}
      \vspace{2ex}
      
      { 
      \noindent Jan Březina, Stanislav Sysala, Jan Stebel,\\[0.6ex]
      Michal Béreš, Simona Bérešová, Pavel Exner,\\[0.6ex]
      David Horák, Jakub Kružík, Tomáš Luber}
      \vspace{2cm}
    
      \noindent leden 2023
    }
    \vfill
% Bottom of the page
    % datum?
\end{titlepage}


\begin{abstract}
Bezpečnost plánovaného hlubinného úložiště radioaktivních odpadů (HÚ) je zajištěna jednak inženýrskými bariérami (obalový soubor, bentonitové utěsnění) a dále zejména geologickou bariérou – horninovým masivem. Pro konstrukci úložiště jsou v~horninovém masivu raženy chodby, přičemž dochází k~narušení okolní horniny a vzniká tzv. zóna porušení (EDZ). V důsledku 
mechanických změn na síti mikro puklin dochází v~EDZ k~nárůstu hydraulické vodivosti. EDZ v~okolí chodeb tak vytváří preferenční cestu pro transport kontaminace ke 
geologickým poruchám -- hlavním transportním cestám k~povrchu. 
Metodika předkládá návrh charakterizace bezpečnosti jednotlivých úložných pozic \sy{z pohledu transportu radioaktivních látek v závislosti na EDZ}. Za tímto účelem je definován indikátor bezpečnosti -- numerická veličina odvozená od transportního modelu a jsou popsány stochastické metody umožňující zachytit vliv nejistot v parametrech modelu. Stochastický model parametrů EDZ je určen pomocí nepřímého měření s využitím Bayesovské inverze.
\end{abstract}

\vfill

      {\noindent {\bf Poděkování:} Výsledek byl realizován s podporou projektu č. TK02010118\\ spolufinancovaného Technologickou agenturou ČR v rámci Programu Théta.
        {\bfseries }\par}

\newpage
\setcounter{tocdepth}{2} % pouze sekce a podsekce 1. úrovně
\tableofcontents


\newpage

\printglossary[type=\acronymtype, nonumberlist, title={Seznam zkratek}]
\glsaddallunused[\acronymtype]
%\glsaddall
\newpage

\section{Úvod}
Nutnou podmínkou povolovacího procesu pro hlubinné úložiště radioaktivních odpadů (HÚ) je prokázání jeho dlouhodobé bezpečnosti. Pro Českou republiku je uvažován koncept HÚ jehož dlouhodobá bezpečnost (milión let) je zajištěna zejména geologickou bariérou a v bližším horizontu (desítky tisíc let) je zajištěna úložným obalovým souborem ÚOS, zejména ocelovými kontejnery a bentonitovým těsněním. Geologická bariéra
zajišťuje nejen zdržení úniku kontaminace na povrch, ale také dostatečné naředění kontaminace. Posouzení bezpečnosti uvedných bariér je realizováno pomocí výpočetních modelů dílčích transportních a difúzních procesů. HÚ pro Českou republiku bude situováno v křehkých horninách, pro které je pórový prostor tvořen prakticky výhradně systémem puklin v široké škále velikostí. 

Existující studie bezpečnosti zahrnují obvykle model uvolňování radionuklidů z matrice uloženého materiálu, dále model porušení kontejnerů a model difúze skrze bentonitovou ucpávku. Následuje model šíření puklinami geosféry a model šíření v~biosféře. Vzhledem ke komplexnosti této problematiky je zde jisté nebezpečí velmi detailního popisu jevů, které je jednoduché popsat (např. inventář radionuklidů, model biosféry), a zanedbání procesů méně známých a složitých. Takovými bezesporu jsou procesy v~geologické bariéře a to zejména kvůli nepřímému pozorování těchto procesů a nejistotám v~popisu geologických struktur. 

V této metodice se věnujeme méně studované části geologické bariéry a sice okolí samotného úložiště, kde dochází k~interakci přirozeného puklinového systému neporušené horniny s modifikovaným systémem puklin v zóně porušení (EDZ -- Excavation Damage Zone) v~okolí přístupových a úložných chodeb. Pro charakterizaci bezpečnosti EDZ v různých částech HÚ je představen model transportu kontaminace od rozhraní bentonit -- hornina, skrze EDZ podél inženýrských struktur až do významnějších hydrogeologických poruch, které mohou sloužit jako transportní cesty k~povrchu. 

Vzhledem k~těmto lokalizovaným transportním cestám a vlivem převažujícího směru proudění podzemní vody lze předpokládat sníženou bezpečnost některých částí úložných chodeb. Uvedený transportní model je proto použit k~definici odvozené numerické veličiny -- indikátoru bezpečnosti jednotlivých pozic ÚOS. Zásadním přínosem metodiky je orientace na použití stochastických metod v~kombinaci s detailními modely transportu a modely vzniku EDZ. To umožňuje predikci indikátoru bezpečnosti se zahrnutím nejrůznějších nejistot v~parametrech transportního modelu. 

Pro praktické použití výpočetních modelů je nutné z vhodných měření stanovit parametry těchto modelů. Jen některé parametry je možné měřit přímo ostatní jsou určovány pomocí inverzních metod. Při použití stochastických výpočtů však potřebujeme znát pravděpodobnostní rozdělení vstupních parametrů. Metodika se proto věnuje i aplikaci Bayesovské inverze pro stanovení rozdělení vstupních parametrů.

Ve zbylé části úvodní kapitoly stručně představíme problematiku zón poškození způsobených ražbou, strukturu hlubinného úložiště radioaktivních odpadů, definici indikátoru bezpečnosti a koncept navrhované metodiky, jehož jednotlivé části budou podrobněji popsány v~následujících kapitolách.



%Stručná rešerše.

%Přehled metodiky. Indikátor, MLMC, Transportní model, Homogenizace, mechanický model, 
%Bayes. 

% Důležitost EDZ je všeobecně akceptována, viz Pusch (2008), a existuje řada
% experimentálních i teoretických prací, které se snaží o~charakterizaci EDZ z hle-
% diska mechanického poškození a návazných hydraulických vlastností, Lanyon
% (2011), Chandler et al. (2002), Vavro (2016).
% Cílem projektu Endorse je vytvoření metodiky pro predikci odvozených ve-
% ličin charakterizujících bezpečnost EDZ, tzv. indikátorů bezpečnostních funkcí
% (dále jen indikátory bezpečnosti). Pro tento účel budeme pojmem EDZ souhrnně
% označovat širší oblast, na které může dojít ke změnám hydraulické vodivosti, viz.
% členění v sekci 3. Výpočet indikátorů bezpečnosti bude založen na modelu trans-
% portu skrze EDZ, jehož parametry budou určeny pomocí modelu vzniku EDZ
% na jemnější prostorové škále. Pro stanovení volných parametrů modelů budou
% na základě rešerše doporučeny vhodné experimentální a měřící metody. Budou
% studovány stochastické výpočetní metody pro aproximaci rozdělení pravděpo-
% dobnosti indikátorů vzhledem k~principiálním nejistotám v popisu horninového
% prostředí. Vyvinuté matematické modely a výpočetní software bude možné vy-
% užít k~obecnějšímu modelování procesů v EDZ, odhadu citlivosti na parametry
% a následně tak i optimalizaci měření nebo monitoringu.
% Tento pracovní dokument popisuje předběžný návrh metodiky a přesněji
% definuje indikátory bezpečnosti. Dále jsou vytipovány hlavní otevřené problémy,
% na které se musí projekt zaměřit, a je shrnuta základní rešerše zdrojů k~dané
% problematice.
% V dokumentu jsou definovány indikátory bezpečnosti EDZ včetně modelů,
% na kterých má být založen jejich výpočet. Pro navržené modely jsou popsány
% 1
% klíčové vstupní parametry a je představena základní rešerše existujících měřících
% metod pro jejich stanovení. Ve spolupráci s aplikačním garantem pak proběhne
% upřesnění experimentů pro validaci metody predikce indikátorů, která bude vý-
% stupem projektu.

\subsection{Zóny poškození způsobené ražbou (EDZ)}

V okolí stěn ražených hlubinných tunelů nebo vrtů v~krystalických horninách jsou pozorovány zvýšené koncentrace napětí, tvorba a spojování trhlin různých rozměrů nebo výlomy ve tvaru V-zářezů. Studium vzniklých zón poškození je nezbytné pro zajištění stability různých důlních děl obecně, nejen pro HÚ. \cite{perras_predicting_2016} navrhli následující členění zóny poškození, viz \figref{fig:edz_zones}:
\begin{figure}
	\centering
	\includegraphics[width=\textwidth]{graphics/EDZ_structure.pdf}
	\caption{Označení zón s kvalitativně odlišnými změnami v~důsledku vzniku díla.}
	\label{fig:edz_zones}
\end{figure}
\begin{description}
	\item[CDZ -- Construction Damage Zone] (\textit{zóna způsobena technologii ražby}).\\ 
    V~případě použití razícího stroje (TBM – Tunnel Boring Machine) je rozsah této zóny zanedbatelný, v~jednotkách centimetrů. V~případě ražby \sy{trhacími prácemi} je rozsah CDZ významný, až 1 m. V~rámci navrhované metodiky budeme způsob ražby zohledňovat pouze volbou parametrů Hoek-Brownova modelu, viz sekce \ref{sec_HB}.
	\item[HDZ -- High Damage Zone] (\textit{zóna s vysokým stupněm poškození}).\\
    Dochází zde k~tvorbě a~spojování puklin na makroskopické úrovni, k~výlomům nebo drcení horninového masivu, přičemž tyto jevy nesouvisí s technologii ražby. 
	\item[EDZ -- Excavation Damage Zone] (\textit{zóna s nevratným poškozením}).\\
    V~této zóně dochází k~nevratným deformačním procesům, například k~tvor\-bě a~spojování trhlin v~menším než makroskopickém měřítku. Podle velikosti trhlin můžeme EDZ dále dělit na vnější zónu poškození (ODZ) a~vnitřní zónu poškození (IDZ), viz \figref{fig:edz_zones}.
	\item[EIZ -- Excavation Influence Zone] (\textit{zóna ovlivněna ražbou}).\\
    Pozorujeme zde zvýšené koncentrace napětí, při kterých ale nedochází k~trvalým deformacím.
\end{description}

V dalším nebudeme pro jednoduchost rozlišovat jednotlivé typy zón a~zkratku EDZ budeme používat pro sjednocení zón HDZ, EDZ a~EIZ. Mezi hlavní faktory, které způsobují tvorbu EDZ, patří počáteční napěťový stav v~horninovém masivu, konkrétně jeho velikost, poměry mezi hlavními napětími a~orientace hlavních směrů napětí vůči ose tunelu. Dále je potřebné zmínit tvar tunelu a~jeho rozměry, způsob ražby, mechanické vlastnosti horninového masivu nebo jeho geologickou strukturu. EDZ se může rozšiřovat i po ražbě, a~to v~důsledku okolní ražby nebo teplotního zahřívání, viz \cite{andersson_aspo_2009}. Naopak k~jejich stabilizaci může například přispět bobtnání bentonitové bariéry HÚ, viz \cite{andersson_aspo_2009}. K dalším změnám pak dochází v~dlouhodobém horizontu po uzavření úložiště a~zatopení díla, viz např. \cite{pusch_can_2012}.

Pro bezpečnost HÚ je kromě mechanické stability důlních děl nezbytné zkoumat také hydro-mechanické procesy v~EDZ, neboť zde dochází ke změnám permeability, porozity či rozevření puklin. Na rozdíl od prostoru chodeb, které mohou být utěsněny různými inženýrskými postupy (bentonit, beton), \sy{je obtížné} EDZ v~rozsahu celého úložiště efektivně utěsnit a~představuje tak hlavní transportní cestu mezi zdrojem kontaminace a~případnými geologickými poruchami protínajícími chodby díla. Námi navržená metodika bude hlavně zaměřena na numerickou analýzu hydro-mechanických procesů v~EDZ a~souvisejí transportní procesy v~celém masivu.

Za účelem predikce EDZ a~analýzy sdružených procesů v~těchto zónách byly a~jsou prováděny různé in-situ experimenty v~podzemních laboratořích. Jako příklad uveďme Äspö pillar stability experiment provedený ve Švédsku, viz \cite{andersson_aspo_2009}, nebo Tunnel Sealing Experiment (TSX) v~Kanadě, viz \cite{Rutqvist2009Modeling}. Připomeňme, že i v~Česku máme takovou laboratoř v~lokalitě Bukov, ve které jsou mimo jiné plánovány experimenty související s~EDZ.

\subsection{Struktura úložiště v ČR}

Pracovní návrh struktury úložiště dle dokumentu SÚRAO\footnote{dokument SÚRAO: \uv{Profily jednotlivých chodeb v HU.pdf}, patrně vychází ze zprávy "Špinka a kol., SÚRAO, 2017", která není veřejně dohledatelná.} představuje základní sché\-ma úložiště, \figref{fig:schema_huro}. Radioaktivní odpad je ukládán v~jedné horizontální rovině v~systému paralelních chodeb kolmých na páteřní chodbu. Délka úložných chodeb/vrtů je několik set metrů. 

V případě horizontálního ukládání jsou jednotlivé ÚOS ukládány roboticky do horizontálních vrtů ražených pomocí razícího stroje (TBM). Vrty jsou navrženy o~průměru 2.2~m s osovou roztečí 20~m. Rozteč ÚOS není přesně stanovena, ale lze předpokládat rozteč minimálně 7~m. Horizontální ukládání má z hlediska vlivu EDZ na bezpečnost výhodu šetrné ražby a~menšího průměru ražby, což obojí přispívá k~menšímu rozsahu EDZ i k~menšímu ovlivnění transportních parametrů. V~případě výrazně rozdílného vertikálního a~horizontálního napětí může být nevýhodný kruhový profil oproti ideálnímu eliptickému profilu.

V případě vertikálního ukládání jsou z páteřní chodby kolmo vedeny ukládací chodby ražené cyklickou metodou drill and blast. V nich jsou pak vrtány svislé ukládací studny. Rozteč úložných chodeb je 18~m, šířka 4~m, výška 6.7~m, předpokládá se přibližně kruhový profil stropu. Ukládací studny o~průměru 1.8~m a~hloubce 6.5~m nemají v~návrhu definovanou osovou rozteč. Vertikální ukládání je patrně technologicky jednodušší, avšak i za předpokladu použití šetrné ražby bude rozsah konstrukční zóny porušení výrazně větší oproti kontinuální ražbě pomocí TBM. Vzhledem k~většímu celkovému průřezu a~nevhodnému profilu ukládací chodby lze předpokládat výrazně větší rozsah EDZ a~vyšší vliv na transportní parametry horniny. 

Zejména kvůli velmi významnému rozsahu CDZ v~případě cyklické ražby a~též kvůli složitější geometrii dále uvažujeme pouze případ horizontálního ukládání. Základní koncept metodiky je však možné aplikovat i pro případ vertikálního ukládání. 

%\todo{bylo by vhodné provést plastické kalkulace pro parametry z TSX a skutečné profily %chodeb, porovnání poškození}
 
\begin{figure}
    \centering
    \includegraphics[width=\textwidth]{graphics/schema_HURO.pdf}
    \caption{Schéma struktury hlubinného úložiště radioaktivního odpadu.}
    \label{fig:schema_huro}
\end{figure}



\subsection{Návrh indikátoru bezpečnosti HÚ založený na transportních procesech}
\label{sec:indikator}

%\sy{[Standa: Navrhl jsem v této sekci několik drobných upřesnění, viz červeně vyznačené texty. Upozornil bych ještě, že v této sekci mluvíme o jedné rovině $\mathcal H$, zatímco v sekci 2.5 o %dvou rovinách.]}
%JB: zapracováno

Vzhledem k~velké celkové ploše HÚ je nutné uvažovat protnutí hydraulicky významných geologických poruch, viz fialová puklina na Obr. \ref{fig:transport_koncept}. V~okolí těchto poruch se předpokládá vynechání úložných pozic ÚOS pro zachování bezpečné doby zdržení a~dostatečného ředění kontaminace v~geologické bariéře. Pro ohodnocení relativní bezpečnosti úložných pozic jednotlivých ÚOS z~hlediska vlivu EDZ byl v~rámci projektu navržen indikátor bezpečnosti úložné pozice (IBUP). Pro uvažovaný únik z dané úložné pozice $p$ (zelený válec na Obr. \ref{fig:transport_koncept}) je uvažován model transportu kontaminace (zelená trasa) skrze EDZ a~okolní horninu do identifikované geologické poruchy a~podél ní až do předepsané vertikální vzdálenosti $H_z$ od horizontální roviny úložiště (horní a~spodní okraj oblasti). 

\begin{figure}
    \centering
    \includegraphics[width=0.75\textwidth]{graphics/EDZ_transport_concept.pdf}
    \caption{Konceptuální model transportu z jedné úložné pozice v~blízkosti HÚ.}
    \label{fig:transport_koncept}
\end{figure}

Pro pevně zvolenou úložnou pozici $p$ je výsledkem transportního modelu pole koncentrací $c_p(t,\vc x)$ závislé na čase a~prostoru. 
Jako indikátor bezpečnosti $I(p)$ je uvažována maximální koncentrace na množině $\mathcal H$ přes všechny kladné časy.
Množina $\mathcal H$ sestává ze dvou horizontálních rovin ve vertikální vzdálenosti $H_z$ od roviny úložiště obsahujícího pozici $p$. 
Položíme-li počátek souřadné soustavy do roviny úložiště, můžeme zapsat definici indikátoru $I(p)$ následovně:
\begin{equation}
    \label{eq:def_indicator}
    I(p) = \max_{t\in (0,\infty)} \max_{\vc x\in \mathcal H} c_p(t, \vc x),\qquad \mathcal H = \{ \vc x : |x_z| = H_z\}, 
\end{equation}
Přesná volba vzdálenosti $H_z$ není pro relativní porovnání bezpečnosti úložných pozic podstatná. 
Výpočty uvedené v~kapitole \ref{sec:transport_numeric_} uvažují $H_z = 30\,\unit{m}$, tedy řádově odpovídající rozteči úložných horizontálních vrtů.

Pro omezení vlivu lokálních odchylek numerického řešení je prakticky
použita aproximace maxima přes množinu $\mathcal H$, která bude diskutována dále při popisu numerického schématu pro transport, viz kapitola \ref{sec:transport_edz}.

Indikátor bezpečnosti bude navíc v~rámci navržené metodiky uvažován jako náhodná veličina ovlivněná řadou nejistot v~parametrech použitých modelů. Hlavní nejistoty jsou naznačeny na Obr. \ref{fig:transport_koncept} oranžově a~patří mezi ně zvýšená vodivost EDZ a~náhodná síť puklin v~hornině. Pomocí víceúrovňové metody Monte Carlo jsou pro náhodnou veličinu indikátoru odhadnuty: střední hodnota a~ rozptyl. Dále je pomocí metody maximální entropie aproximována hustota pravděpodobnosti indikátoru. Podrobný popis obou metod je uveden v~kapitole \ref{sec:mlmc}.

\sy{
Navržený indikátor bezpečnosti lze použít pro relativní srovnání bezpečnosti úložných pozic, avšak z hlediska jeho interpretace by bylo lépe definovat ho ve smyslu výsledné dávky ozáření. K tomu by však bylo nutné nahradit vhodnou aproximací navazující transport celou geologickou bariérou, dále šíření biosférou a uvažovat inventář radionuklidů s dlouhým poločasem rozpadu. V rámci projektu byla úspěšně testována aproximace 1D transportu v homogenním prostředí se sorpcemi a rozpady. Z dřívějších studií je známo, že komplexní 3D transportní model lze poměrně dobře aproximovat tímto zjednodušeným 1D modelem jehož parametry je však třeba kalibrovat pomocí plného 3D modelu. Provedení těchto kalibrací a sestavení dalších nutných aproximací přesahovalo časový rámec projektu nehledě k problematice vstupních dat. 
}

Ačkoliv je vlastní indikátor bezpečnosti navržen pro relativní porovnání bezpečnosti jednotlivých pozic $p$, použité výpočetní modely lze po zahrnutí sorpce a~rozpadů případně dále využít jako vstup do modelů transportu v~geologické bariéře.

Pro puklinové prostředí je charakteristická vysoká heterogenita hydraulické vodivosti, která je na makroskopických puklinách i o~několik řádů vyšší nežli vodivost kontinua, zahrnující pukliny menších škál. Pro zachycení tohoto kontrastu je využita kombinace kontinuálního popisu a~sítě diskrétních puklin, kterou budeme v~dalším označovat zkratkou DFM (discrete fracture -- matrix).

\subsection{Idea navrhované metodiky a členění dokumentu}

% \sy{[Standa: Pokusil jsem se tuto sekci rozepsat podrobněji, aby bylo patrnější, jaké jsou hlavní myšlenky navržené metodiky. Také jsem navrhl změnu v~organizaci textu, viz kapitoly 5, 6 a 7. Prosím o~revizi červeně vyznačeného textu a o~odsouhlašení navrženého členění textu.]}
%JB: OK - drobné úpravy

Pro posouzení vlivu EDZ na výše uvedený bezpečnostní indikátor navrhujeme kombinaci matematických modelů a~metod, které zahrnují modely transportu a~poroelasticity v~horninovém prostředí bez puklin i s puklinovou síti, homogenizační techniky, stochastické metody, inverzní analýzu a~další pokročilé metody. V~následujících odstavcích popíšeme hlavní myšlenky navrhovaného konceptu metodiky a~členění tohoto dokumentu.

Kapitola \ref{sec:transport_edz} obsahuje popis modelu transportu v~blízkém okolí HÚ se zahrnutím náhodných a~předepsaných puklin. Tento model  je klíčový pro praktické vyhodnocení bezpečnostního indikátoru dle výše uvedené definice. Dále upřesníme, jaké parametry tohoto modelu mohou ovlivňovat EDZ. Nakonec představíme víceškálový transportní model obsahující homogenizační techniky, díky kterým můžeme řešit rozsáhlejší úlohy.


V kapitole \ref{sec:poroelasticity} představíme modely poroelasticity bez puklinové sítě i s puklinovou síti a~příslušné numerické metody. Jedná se o~základní modely, které budou používány pro studium hydro-mechanických procesů v~EDZ. Tyto modely lze použít k~numerické simulaci plánovaných nebo již provedených in-situ experimentů. Na základě experimentálních měření (například pórových tlaků) můžeme  upřesnit materiálové parametry poroelastických modelů, a~to pomocí inverzní analýzy. Takto upřesněné parametry nakonec využijeme i v~transportním modelu, díky čemuž zahrneme vliv EDZ na hodnotu bezpečnostního indikátoru.

Nedílnou součástí navrhované metodiky je výzkum nejistot pomocí stochastických metod, které představíme v~kapitole \ref{sec:nejistoty}. Nejprve popíšeme víceúrovňovou metodu Monte Carlo, kterou používáme k~propagaci nejistot skrze transportní model a~výpočetně efektivnímu odhadu pravděpodobnostního rozdělení indikátoru bezpečnosti. Následně popíšeme Bayesovskou inverzní analýzu, která umožňuje určit pravděpodobnostní rozdělení neznámých vstupních parametrů na základě naměřených dat zatížených nejistotami (šumem). Bayesovskou inverzní analýzu budeme především aplikovat na úlohu poroelasticity.

{\sloppy
Praktická aplikace navržené metodiky je podpořena vývojem specializovaného softwarového nástroje Endorse, který ji umožňuje realizovat. V~kapitole \ref{sec:SW} proto stručně popíšeme schéma tohoto software a~jeho hlavní komponenty.
}

V kapitole \ref{sec:transport_numeric_} ilustrujeme navržený koncept metodiky na numerickém příkladu, přičemž provádíme predikci indikátoru bezpečnosti pomocí víceúrovňové metody Monte Carlo a~transportního modelu. V kapitole \ref{sec:TSX} popíšeme část TSX experimentu realizovaného v~Kanadě (Underground Research Laboratory), která je relevantní navržené metodice. Na základě naměřených dat pak použijeme Bayesovskou inverzní analýzu k~upřesnění parametrů poroelastických modelů popisujících změny hydraulické vodivosti v EDZ.
V~závěrečné kapitole \ref{sec:zaver} shrneme hlavní principy představené metodiky. 

V dodatku \ref{sec:rozsireni} diskutujeme další možnosti rozšíření metodiky, které zahrnují numerickou predikci rozsahu plastických změn v~EDZ, 
užití pokročilých numerických metod řešení poroelastických modelů s puklinami i bez nich a~definování hydraulické vodivosti na poli posunutí pomocí homogenizačních technik.
Vzhledem k~absenci komplexních experimentálních dat nebylo možné metodiku plně validovat. V dodatku \ref{sec:Bukov} proto přikládáme návrh experimentu 
za účelem ověření prezentované metodiky.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model transportu \sy{radionuklidů} v okolí úložiště}
\label{sec:transport_edz}


Pro výpočet indikátoru bezpečnosti je uvažován transport konzervativního stopovače skrze okolní EDZ a~dále skrze deterministické (známé) a~stochastické (neznámé) pukliny v~okolí HÚ. 
Zdrojem koncentrace je tok koncentrace z~bentonitové obálky jednoho úložného místa předepsaný pomocí koncentrace uvnitř obálky a~difúzního odporu obálky. 
Použitý model proudění a~transportu v~kontinuálním popisu je popsán v~sekci \ref{sec:transport_rovnice}. Následně v~části \ref{sec:vek_pocet_dfm} zavedeme značení pro DFM popis a~v~části \ref{sec:transport_DFN} popíšeme DFM model proudění a~transportu v~puklinovém prostředí. V~sekci \ref{sec:transport_num} je stručně popsána použitá diskretizace transportního modelu použitá v~simulátoru Flow123d 
vyvíjeném na TU Liberec, viz \cite{Brezina2021Flow123d}. V~části \ref{sec:transport_indikator} aproximujeme výše uvedenou definici indikátoru bezpečnosti $I(p)$ tak, aby byla robustní vůči navrženému numerickému řešení transportního problému. Poslední sekce \ref{sec:transport_multiscale} diskutuje numerickou homogenizaci použitou pro tvorbu hrubších aproximací modelu nutných k~propagaci nejistot pomocí víceúrovňové metody Monte Carlo (MLMC). Samotná metoda MLMC bude popsaná později v~sekci \ref{sec:mlmc}.

Jelikož indikátor bezpečnosti má sloužit pro relativní porovnání jednotlivých úložných míst, není pro jednoduchost uvažována sorpce ani radioaktivní rozpad.
Implementace modelu transportu v~simulátoru Flow123d ovšem oba jevy zahrnuje, včetně simulace transportu a~rozpadové řady více radionuklidů. Použití reálného inventáře radionuklidů je proto možné, pokud budou dostupné příslušné parametry pro všechny izotopy. 

\subsection{Transportní rovnice}
\label{sec:transport_rovnice}
Transport kontaminace je způsoben prouděním podzemní vody, které je pro malé rychlosti popsáno rovnicí kontinuity a~Darcyho zákonem.
Pro dlouhodobě ustálené podmínky lze uvažovat stacionární podobu rovnic proudění ve tvaru:
\begin{equation}
    \label{eq:darcy_matrix}
    \div \vc v = 0, \quad \vc v = -\tn K \grad (p + z),
\end{equation}
kde $\vc v$ $\bunit{m\cdot s^{-1}}$ je Darcyho (makroskopická) rychlost proudění, $\tn K$ $\bunit{m\cdot s^{-1}}$ je tenzor hydraulické vodivosti, $p$ $\bunit{m}$ je tlaková výška a $z$ $\bunit{m}$ je vertikální složka polohového vektoru. Při modelování transportu skrz EDZ hraje důležitou roli zejména tenzor hydraulické vodivosti, který může být značně heterogenní a anizotropní. Vliv EDZ na $\tn K$ bude popsán v~kapitole \ref{sec:poroelasticity} s využitím poroelastického modelu.

Po vyřešení rovnice \eqref{eq:darcy_matrix} s vhodnými okrajovými podmínkami a získáním rychlostního pole  $\vc v$ můžeme časový vývoj koncentrace konzervativního stopovače $c$ $\bunit{kg\cdot m^{-3}}$ popsat rovnicí advekce-difúze:
\begin{equation}
   \label{eq:transport_matrix}    
   \prtl_t (\theta c) + \div(\vc j) = 0,\quad \vc j = \vc v c - \theta \tn D \grad c,
\end{equation}
kde $\theta$ $[-]$ je porozita, $\tn D$ $\bunit{m^2\cdot s^{-1}}$ je tenzor difúze a hydrodynamické disperze, dále jen {\it disperzní tenzor}. Difúze je způsobena chaotickým pohybem molekul roztoku, zatímco disperze je důsledkem nehomogenity rychlostního pole na mikroskopické úrovni, viz sekce \ref{sec:transport_multiscale}.  Pro případ izotropního média je disperzní tenzor dán vztahem:
\begin{equation}
  \label{eq:dispersion_tn}
  \tn D = \tn D_m + \abs{\vc u}\Big(\alpha_T \tn I + (\alpha_L - \alpha_T)\frac{\vc u \otimes \vc u}{\abs{\vc u}^2}\Big),
\end{equation}
kde 
\begin{itemize}
 \item $\tn D_m = d_m \tn I$ je tenzor molekulární difúze, který pro jednoduchost uvažujeme izotropní s~difúzním koeficientem $d_m$ $\bunit{m^2\cdot s^{-1}}$, řádově $10^{-9}$ pro čistou vodu,
% \item $\tau=\theta^{1/3}$ $[-]$ je tortuozita dle %\cite{millington_permeability_1961},
 \item $\alpha_L$ a $\alpha_T$ $\bunit{m}$ jsou koeficienty podélné a příčné disperze,
 \item $\vc u$ je pórová rychlost,  $\vc u = \vc v / \theta$.
\end{itemize}
Zdroj stopovače je předepsán pomocí okrajové podmínky na stěně úložného vrtu:
\begin{equation}
    \label{eq:source_bc}
    \vc j \cdot \vc n = f_S = \chi(\vc x) \alpha_S \big( c_S(t) - c(t, \vc x\big)).
\end{equation}
Zde na levé straně je difúzní tok ve směru normály $\vc n$ přes hranici, 
jelikož při výpočtu rychlostního pole je hranice uvažována jako nepropustná a rychlostní pole je kolmé k~normále. Plošná hustota difúzního toku $f_c$
 $\bunit{kg\cdot s^{-1}\cdot m^{-2}}$ na pravé straně je dána charakteristickou funkcí $\chi(\vc x)$ zdroje kontaminace, difúzním parametrem $\alpha_B$ $\bunit{s^{-1}}$
 a referenční koncentrací $c_S(t)$ $\bunit{kg\cdot m^{-3}}$. Difúzní parametr zahrnuje difúzi skrze ÚOS, tedy skrze porušený kontejner a skrze bentonitovou vrstvu, viz \cite{nordman_radionuclide_2022}, sekce 2.2. Referenční koncentrace $c_S(t)$ reprezentuje koncentraci uvnitř kontejneru. Kvůli nedostatku zdrojů uvažujeme exponenciální rozpad $c_S(t) = c_0 \exp(-\alpha_S t)$, coby řešení jednoduché lineární diferenciální rovnice.

Vývoj koncentrace konkrétních radionuklidů je silně ovlivněn též radioaktivním rozpadem a sorpcí. Pro zjednodušení výpočtů a zejména interpretace výsledků není vliv těchto jevů v~metodice uvažován. Simulátor Flow123d implementuje oba jevy v~plné obecnosti. Níže ukážeme, že použité zjednodušení je pro relativní porovnání úložných pozic přijatelné v~případě homogenní porozity $\theta$ a sorpce.

\subsection{Vektorový počet pro DFM}
\label{sec:vek_pocet_dfm}

V této části zavedeme vektorový počet pro tzv. DFM (discrete fracture-matrix) popis puklinového prostředí, kombinující ekvivalentní kontinuum a síť diskrétních puklin.  Označíme $\Omega_m \subset \Real^3$ oblast matrice a  $\Omega_f=\Cup_{i=1}^{n_{fr}}\mathcal F_i$ oblast puklin, přičemž předpokládáme, že jednotlivé pukliny $\mathcal F_i$ mají polygonální tvar.
%$\mathcal F_i: D_i \to \Real^3$ jsou parametricky zapsané polygony. 
% \sy{[Standa: Pro účely tohoto reportu myslím není potřeba psát, že pukliny jsou zapisovány parametricky. Dále zvažte, jestli nezavést značení $\Omega=\Omega_m\cup\Omega_f$, přičemž symbol $\Omega$ bychom používali pro definici kontinuálních modelů, např. v~sekci 2.1, což momentálně chybí.]}
% \jb{Je dobré si uvědomit, že veličiny na puklinách jsou funkce pouze dvou parametrických souřadnic.}
Pro jednoduchost nebudeme dále diskutovat průsečíky puklin, nicméně 
použitá implementace v~simulátoru Flow123d s nimi počítá. 
Pro danou skalární, vektorovou, nebo tenzorovou veličinu $p=(p_m, p_f)$ označíme $p_m$ funkci veličiny na $\Omega_m$ a $p_f$ funkci definující veličinu na $\Omega_f$ (fakticky je funkce definována na parametrických definičních oborech $D_i$). 
Pro každý bod $\vc x \in \Omega_f$ na jednotlivé puklině označíme $\vnu(\vc x)$ normálový vektor definující též její orientaci. Dále budeme rozlišovat kladnou a zápornou stranu plochy $\vnu^+ = \vnu$ a~$\vnu^-=-\vnu$. Veličina $p$ může být v~normálovém směru k~puklině nespojitá, proto pro $p_m$ definujeme její hodnotu na stranách pukliny: $p_m^+$, $p_m^-$. Pole tloušťky pukliny označíme $\delta(\vc x)$.


Pro zjednodušení zápisu rovnic na puklinách nyní zavedeme semi-diskrétní diferenciální operátory. Pro skalární veličinu $p$ je $p_f$ funkcí dvou proměnných a~její tečný gradient $\grad_t p_f$ je pak tečným vektorem k~$\Omega_f$. Semi-diskrétní gradient pak aproximuje normálovou složku pomocí diference stop:
$$
\fgrad p = \grad_t p_f + \grad_\nu p,\qquad 
\grad_\nu p = \frac{1}{2}\Big(\Delta^+ p\vnu^+  + \Delta^- p\vnu^-\Big),
$$
přičemž normálová složka je aproximována součtem diferencí přes jednotlivé strany pukliny:
$$
\qquad \Delta^\pm p = \frac{2}{\delta}\Big(p_m^\pm - p_f\Big).
$$
Dále definujeme semi-diskrétní gradient na kladné a záporné straně plochy pukliny:
$$ \fgrad^\pm p = \grad_t p_f + \Delta^\pm p\vnu^\pm. $$
Tato definice umožňuje formální rozšíření notace pro průsečíky puklin.
% \todo{Vypadá to, že rovnice sice jsou zapsat těmito posčítanými operátory, ale pro zápis Neumanovy okrajové podmínky, která dělá vazbu stejně potřebujeme zavést $\fgrad^\pm$.}

Stejným způsobem zavedeme gradient pro vektorové pole $\vc v=(\vc v_m, \vc v_f)$:
$$
\fgrad \vc v = \grad_t \vc v_f + \grad_\nu \vc v,\qquad 
\grad_\nu \vc v = \frac{1}{2}\Big(\Delta^+\vc v \otimes\vnu^+  + \Delta^- \vc v\otimes\vnu^-\Big),
$$
Analogicky pak zavedeme semi-diskrétní divergenci vektorového pole:
$$
\fdiv \vc v = \div_t \vc v_f + \div_\nu \vc v,\qquad
\div_\nu \vc v = \frac{1}{2}\Big(\Delta^+ \vc v \cdot \vnu^+  + \Delta^- \vc v \cdot \vnu^-\Big)
$$
Zde tečná divergence $\div_t \vc v_f$ je divergencí tečné složky pole $\vc v_f$. Podobně pak zavedeme též divergenci pro tenzorovou veličinu $\vc\sigma = (\vc\sigma_m, \vc\sigma_f)$:
$$
\fdiv \vc\sigma = \div_t \vc\sigma_f + \div_\nu \vc\sigma,\qquad
\div_\nu \vc\sigma = \frac{1}{2}\Big(\Delta^+ \vc\sigma  \vnu^+  + \Delta^- \vc\sigma  \vnu^-\Big).
$$



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Transportní rovnice pro DFM popis}
\label{sec:transport_DFN}

Pro DFM popis je rovnice \eqref{eq:darcy_matrix} uvažována pouze na 
oblasti kontinua $\Omega_m$. Na puklinové síti $\Omega_f$ je proudění popsáno rovnicí:
\begin{equation}
    \label{eq:darcy_fracture}
    \fdiv (\delta \vc v_f) = 0, \quad \vc v_f = -\tn K_f \fgrad (p + z).
\end{equation}
Zde $\tn K_f(\vc x)$ je 3D tenzor vodivosti popisující materiál puklinové sítě. Vazba mezi oblastmi je dána bilancí objemového toku na rozhraních mezi matricí a puklinou, což představuje Neumanovu okrajovou podmínku pro rovnici \eqref{eq:darcy_matrix} na $\Omega_m$. Je třeba uvažovat samostatnou podmínku pro každou stranu pukliny:
\begin{equation}
    \vnu^{\pm} \cdot \vc v_m^\pm = \vnu^{\pm} \cdot \vc v_f^\pm =  
    -\vnu^\pm \cdot (\tn K_f \fgrad^\pm (p + z))
    \quad \text{na }\Omega_f.
\end{equation}
Pro model proudění uvažujeme Dirichletovu okrajovou podmínku pro piezometrickou výšku $P = p + z$, která je interpolována z~regionálního modelu proudění na modelové lokalitě. 

Analogicky, pro DFM popis transportu uvažujeme 
rovnici \eqref{eq:transport_matrix} pouze pro matrici $\Omega_m$,
zatímco na puklinách $\Omega_f$ má transportní rovnice tvar:

\begin{equation}
   \label{eq:transport_fr}    
   \prtl_t (\delta\theta c_f) + \fdiv \,\vc j = 0, \quad  
   \vc j_f = \vc v_f c_f - \delta \theta \tn D \fgrad c.
\end{equation}
Vazba mezi oblastmi je podobně jako pro proudění dána bilancí hmotnostního toku stopovače:
\begin{equation}
    \vnu^{\pm} \cdot \vc j_m^\pm  = \vnu^{\pm} \cdot \vc j_f^\pm  \quad \text{na }\Omega_f.
\end{equation}
Na vnějším okraji modelu je uvažován nulový hmotnostní tok na vtokové části hranice a nulový difúzní tok na výtokové části hranice.


\subsection{Numerické řešení transportních úloh}
\label{sec:transport_num}

Numerické řešení transportních úloh je realizován pomocí software Flow123d. Konkrétně rovnice proudění \eqref{eq:darcy_matrix} a \eqref{eq:darcy_fracture} jsou řešeny pomocí smíšené hybridní metody s využitím Raviartových-Thomasových prvků nultého řádu pro diskretizaci rychlostního pole a po částech konstantních prvků pro diskretizaci tlaku. Kontinuita toků je vynucena Lagrangeovými multiplikátory na stěnách prvků. Výpočetní oblast je rozdělena na simplexové elementy tak, aby puklinové 2D elementy byly totožné se stěnami přilehlých 3D elementů. Diskretizovaný problém pak má algebraickou reprezentaci
% 
\eq{\label{eq:flow_alg_mh} 
\begin{bmatrix}
\tn F^{\vc v\vc v} & \tn B^{\vc v p} & \tn B^{\vc v\lambda}\\
(\tn B^{\vc v p})^\top & \tn O & \tn O\\
(\tn B^{\vc v\lambda})^\top & \tn O & \tn O
\end{bmatrix}
%
\begin{bmatrix}
\uu^{\vc v}\\ 
\uu^p\\
\uu^\lambda
\end{bmatrix} = 
%
\begin{bmatrix}
\vc 0\\
\vc f^p\\
\vc 0
\end{bmatrix}, 
}
% 
kde $\uu^{\vc v}$, $\uu^p$ a $\uu^\lambda$ jsou stupně volnosti pro pole rychlosti, tlakové výšky v~těžištích elementů výpočetní sítě a pro Lagrangeovy multiplikátory. % \sy{[Pozn. od EH: doplnila bych vysvětlující význam pravých stran soustavy.]}
Vektor $\vc f^p$ je diskretizovaný zdrojový člen v bilanci hmoty.
Matice $\tn F^{\vc v\vc v}$ a $\tn B^{\vc vp}$ jsou sestaveny lokálně po elementech, což umožňuje efektivní výpočet Schurova doplňku. Soustavu lineárních rovnic \eqref{eq:flow_alg_mh} tak lze redukovat na soustavu pro $\uu^\lambda$ se symetrickou pozitivně definitní maticí, která je řešena předpodmíněnou metodou sdružených gradientů.

Transportní rovnice \eqref{eq:transport_matrix} a \eqref{eq:transport_fr} jsou řešeny pomocí nespojité Galerkinovy metody (DGM). Obvykle používáme prvky prvního nebo druhého řádu.
Pro vyšší robustnost vzhledem k~poměru advekce a difúze (disperze) je zvolena nesymetrická varianta DGM, která vede na algebraickou soustavu
$$ \tn T\uu^c = \vc f^c $$
s obecně nesymetrickou maticí. Soustava je řešena metodou bikonjugovaných gradientů s vhodným předpodmíněním. 

\subsection{Aproximace indikátoru bezpečnosti}
\label{sec:transport_indikator}

V sekci \ref{sec:indikator} byl zaveden indikátor bezpečnosti $I(p)$ konkrétní úložné pozice $p$ reprezentující jeden ÚOS, viz rovnice \eqref{eq:def_indicator}. Indikátor určujeme na základě znalostí pole koncentrací $c(t, \vc x)$, které vypočítáme pomocí transportního modelu definovaného výše v této kapitole. Pole koncentrací závisí na pozici zdroje kontaminace dané charakteristickou funkcí $\chi(\vc x)$ v okrajové podmínce \eqref{eq:source_bc}.

Abychom zdůraznili závislost pole koncentrací na úložné pozici $p$, používáme v rovnici \eqref{eq:def_indicator} značení $c_p(t, \vc x)$ namísto $c(t, \vc x)$. Dále připomínáme, že při definici $I(p)$ vyhodnocujeme pole koncentrací ve vzdálenosti $H_z$ od roviny úložiště, tj. na dvou rovnoběžných rovinách, které se nacházejí nad a pod rovinou s horizontálními úložnými vrty, viz \figref{fig:indicator_cuts}. Konkrétně, na těchto dvou rovnoběžných rovinách hledáme maximum $c_p(t, \vc x)$ přes zkoumaný časový interval a prostorovou proměnnou $\vc x$.

\begin{figure}[ht]
    \centering
    \includegraphics[width=\textwidth]{graphics/indicator_cuts_conc_contours_3D.png}
    \caption{Roviny pro vyhodnocení indikátoru bezpečnosti. Na horní rovině je zobrazen řez polem koncentrace s izoliniemi v~čase maximální hodnoty indikátoru.}
    \label{fig:indicator_cuts}
\end{figure}

Pro danou úložnou pozici $p$ je zvolena oblast $\Omega_p$ zahrnující okolí této pozice, na které počítáme pole koncentrací $c_p$. Tato oblast, včetně diskrétních puklin v ní obsažených, je následně vysíťována, viz \figref{fig:indicator_cuts}. Volba velikosti výpočetní oblasti a jemnosti její diskretizace je kompromis mezi přesností výpočtu a výpočetní náročností. Pro orientační výpočet používáme oblast ve tvaru kvádru o~výšce a~šířce 60~m a vodorovné délce 200~m. Geometrie úlohy bude podrobněji popsána v kapitole \ref{sec:transport_numeric_}. Pro zahrnutí vlivu okolních struktur HÚ 
bude vhodné volit větší oblast.

{\sloppy
Vzhledem k~numerickému výpočtu pole koncentrací $c_p$ je vhodné přistoupit k~aproximaci definice indikátoru bezpečnosti, a to ze dvou důvodů. Za prvé, transportní rovnice je řešena pomocí nespojité Galerkinovy metody, která má řadu výhod, avšak výsledná aproximace pole koncentrací $c_p$ nemusí být v~prostoru spojitá. Konkrétně pro výpočet $I(p)$ používáme po částech konstantní aproximaci $c_p$, která vznikne vyhodnocením nespojitého řešení v~těžištích jednotlivých konečných prvků. Tímto výrazně redukujeme vyšší koncentrace $c_p$, které se mohou vyskytnout na hranách konečných prvků v důsledku numerických chyb. 
}

Za druhé, potřebujeme vyhodnotit pole $c_p$ na horizontálních řezech ve vzdálenosti $H_z$ od roviny úložiště. Poloha rovin $H_z$ je volena tak, aby tyto roviny byly blízko hranice oblasti $\Omega_p$, ale stále uvnitř $\Omega_p$. Hodnoty na elementech protnutých rovinami jsou uloženy do pole a z~tohoto pole hodnot je vyhodnocen kvantil $1-q$ na hladině \sy{významnosti} $q$. Pro výpočetní síť okolo 100 tisíc elementů je vhodná hodnota $q=0.001$. Získaná časová sekvence kvantilů je interpolována kubickým splinem, z~něhož je následně vyhodnoceno maximum přes časový interval simulace.






%V sekci \ref{sec:indikator} byl zaveden indikátor bezpečnosti $I(p)$ konkrétní úložné pozice $p$
%jako maximum přes prostor a čas z pole koncentrace $c(t, \vc x)$ na dvojici vodorovných rovin ve vertikální vzdálenosti $H_z$ od roviny úložiště, viz rovnice \eqref{eq:def_indicator}.
%V této sekci diskutujeme praktické vyhodnocení indikátoru na základě numerického transportního modelu. Pro danou úložnou pozici $p$ je zvolena výpočetní oblast $\Omega_p$ zahrnující okolí úložné pozice. Volba velikosti této oblasti a jemnosti její diskretizace je kompromis mezi přesností výpočtu a výpočetní náročností. Pro orientační výpočet postačí krychle o~hraně 60~m, nicméně pro zahrnutí vlivu okolních struktur HÚ je vhodné volit větší oblast. 

%Oblast $\Omega_p$ včetně diskrétních puklin je následně vysíťována a pro zdroj kontaminace na úložné pozici $p$ je pomocí transportního modelu vypočteno pole koncentrací $c_p(t, \vc x)$.  
%Transportní rovnice je řešena pomocí nespojité Galerkinovy metody, která má řadu výhod, avšak
%výsledné pole koncentrací nemusí být v~prostoru spojité. Jako pole $c_p(t, \vc x)$ je proto použita po částech konstantní aproximace, která vznikne vyhodnocením nespojitého řešení v~těžištích elementů. Tato projekce výrazně redukuje výskyt odlehlých hodnot. 

%Dalším krokem je vyhodnocení pole $c_p$ na horizontálních řezech ve vzdálenosti $H_z$ od roviny úložiště, viz \figref{fig:indicator_cuts}. Poloha rovin $H_z$ je volena tak, aby tyto roviny byly blízko hranice, ale stále uvnitř výpočetní oblasti. Hodnoty na elementech protnutých rovinami jsou uloženy do pole a z tohoto pole hodnot je vyhodnocen kvantil na hladině $1-q$. Pro výpočetní síť  okolo 100 tisíc elementů je vhodná hodnota $q=0.001$. Získaná časová sekvence kvantilů je interpolována kubickým splinem, z něhož je následně vyhodnoceno maximum přes časový interval simulace.


\subsection{Víceškálový transportní model}
\label{sec:transport_multiscale}
Pro zachycení prostorové variability hydraulické vodivosti v~EDZ je třeba rozlišení okolo 0.3~m. Naproti tomu chceme zachytit také vliv sousedních chodeb úložiště. Plné rozlišení EDZ na celé chodbě se zdrojem kontaminace i pro chodby sousedící by vyžadovalo desítky milionů elementů a výpočetní časy jedné realizace
v řádu desítek hodin sekvenčního výpočtu. V~této sekci proto popíšeme postup numerické homogenizace, který umožní přenos jemné prostorové variability parametrů transportního modelu na hrubší výpočetní sítě. Tento postup je případně možno aplikovat ve více úrovních a získat tak sekvenci méně přesných, avšak výpočetně laciných aproximací pro víceúrovňovou metodu Monte Carlo popsanou níže v~sekci \ref{sec:mlmc}. 

\def\macroT{\mathcal{T}_H}
\def\microT{\mathcal{T}_h^E}
Cílem homogenizace je určení ekvivalentních hodnot parametrů modelu na množině elementů makroskopické sítě $\macroT$. Stejná výpočetní oblast je rozdělena na více překrývajících se podoblastí tak, že každý makro element $E\in \macroT$ je celý uvnitř alespoň jedné podoblasti $\Omega_E$. Následně je oblast $\Omega_E$ vysíťována elementy jemné sítě $\microT$. 
Pro element $E$ hrubé sítě $\macroT$ nalezneme podmnožinu elementů jemné sítě $\cS_E \subset \microT$, které se s ním protínají. Zde elementy $\cS_E$ zahrnují
i~puklinové elementy. 
Na jednotlivých elementech jemné sítě $e \in \microT$ uvažujeme heterogenní a po elementech konstantní parametry modelu:
porozitu $\theta^e$, tenzor hydraulické vodivosti $\tn K^e$, disperzivity $\alpha_L^e$, $\alpha_T^e$, a difúzní parametr zdrojového členu $\alpha_S^e$. Každý z~těchto parametrů je třeba homogenizovat odlišným způsobem.


{\it Ekvivalentní} hodnota porozity $\theta_E$ na makro elementu $E$ je vypočtena jako vážený průměr:
\begin{equation}
    \theta_E = \avg{\theta}_E = \frac{1}{|\cS_E|} \sum_{e\in \cS_E} |e| \theta^e,\qquad |\cS_E| = \sum_{e\in \cS_E} |e|,
\label{Eq:average} 
\end{equation}

Ekvivalentní tenzor hydraulické vodivosti $\tn K_E$ nalezneme pomocí požadavku na jeho akci vůči vektorovým polím rychlosti a gradientu tlaku. Na jemné síti $\microT$ řešíme rovnice \eqref{eq:darcy_matrix} a
\eqref{eq:darcy_fracture} s Dirichletovou okrajovou podmínkou pro piezometrickou výšku:
$$
    p(\vc x) + z = \vc U \cdot \vc x.
$$
Pro daný vektor $\vc U$ tak dostáváme po částech konstantní pole rychlostí $v^e(\vc U)$ na elementech jemné sítě $e\in \microT$. 
Pro vektorové pole rychlostí $\vc v^e$ na jemné síti určíme příslušné pole tlakových gradientů $\vc g^e$ z~rovnic:
$$
    \vc v^e(\vc U) = \tn K^e \vc g^e(\vc U).
$$
Tento postup je vhodný pro použitou smíšenou formulaci, kde diskrétní pole tlaků nemá přirozeně definované gradienty.

Nyní definujeme $\tn K_E$ vztahem:
\begin{equation}
    \label{eq:homo_cond_tn} 
    \avg{\vc v^e(\vc U)}_E = \tn K_E  \avg{\vc g^e(\vc U)}_E,\quad  
\end{equation}
kde průměry $\avg{\vc v_k}_E$ a $\avg{\vc g_k}_E$ jsou definovány podle \eqref{Eq:average}. 
Rovnici \eqref{eq:homo_cond_tn} uvažujeme pro tři nezávislé vektory $\vc U$ a výslednou přeurčenou soustavu lineárních rovnic pro 6 hodnot
symetrického tenzoru $\tn K_E$ řešíme ve smyslu nejmenších čtverců.

Určení ekvivalentního tenzoru disperze je výrazně komplikovanější. 
Disperze zachycuje vliv fluktuací rychlostního pole na menších škálách a
je tedy ovlivněna heterogenitou a anizotropií tenzoru hydraulické vodivosti. Vztah \eqref{eq:dispersion_tn} platí pouze pro izotropní prostředí v~obecném případě je závislost disperzního tenzoru $\tn D$ rychlostním poli popsána čtyřrozměrným tenzorem (\cite{bear_tensor_1961}, \cite{scheidegger_general_1961}):
$$
  \label{eq:dispersion_4d}
  \tn D_{mn} = \frac{1}{\abs{\vc u}}\alpha_{mnij} u_i u_j
$$
Přesná podoba numerické homogenizace pro tento obecný případ je 
předmětem dalšího výzkumu. Zatím je proto použit pouze izotropní tenzor vodivosti s~ekvivalentními disperzivitami průměrovanými jako v~případě porozity.

Posledním speciálním případem je homogenizace okrajové podmínky. Hrubá síť není schopna zachytit detailní geometrii úložné chodby, proto je nutné okrajovou podmínku \eqref{eq:source_bc} nahradit objemovým zdrojem:
$$
    f_S^E = \alpha_S^E(c_S(t) - c(t, \vc x)).
$$
Každé okrajové stěně $s$ jemné sítě s aktivním zdrojem ($\chi(s) = 1$)
je přiřazen právě jeden element hrubé sítě $E(s)$, se kterým má největší průnik. Následně je hodnota difúzního koeficientu na makro elementu $E$ určena ze vztahu:
$$
    \alpha_S^E = \frac{1}{V_E}\sum_{E(s)=E} \abs{s} \alpha_S.
$$

% \subsection{Vliv sorpce a rozpadů}

% Sorpci i rozpad jednoho radionuklidu lze zahrnout přidáním členu na pravou stranu rovnice \eqref{eq:transport_matrix}:
% \begin{equation}
%    \label{eq:transport_matrix_decay}    
%    \prtl_t (\theta c) + \div(\vc j) = f - R c.
% \end{equation}
% V případě homogenní sorpce je parameter $R$ konstatní. Dosazením do rovnice snadno ověříme, že řešení rovnice \eqref{eq:transport_matrix_decay} má tvar:
% $$
% c(t, \vc x) = \exp(-\frac{R}{\theta}t) \tilde c(t, \vc x),
% $$
% kde $\tilde c(t, \vc x)$ je řešení původní rovnice \eqref{eq:transport_matrix}
% s přeškálovaným zdrojovým členem v okrajové podmínce:
% $$
% \tilde f_S = f_S / \exp(-\frac{R}{\theta} t)
% $$. 

% Jelikož je rovnice lineární bude i  

% Uvažujme konstantní parametr $R$, což je platné pro radioaktivní rozpad, ale nemusí být přijatalné pro silnou hetrogenitu sorpcí. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Model poroelasticity a jeho numerické řešení}
\label{sec:poroelasticity}


Hodnota indikátoru bezpečnosti, definovaného v~sekci \ref{sec:indikator}, je ovlivněna transportní rychlostí $\vc v$, která je přímo úměrná hydraulické vodivosti a gradientu piezometrické výšky. V~okolí úložného místa je kvůli přítomnosti EDZ hydraulická vodivost výrazně vyšší oproti intaktnímu masivu.
Změna vodivosti souvisí se vznikem puklin a jejich otevíráním vlivem mechanického namáhání v~okolí vyraženého díla. 

Abychom mohli odhadnout změnu permeability v~EDZ, je zapotřebí uvažovat nejen porézní proudění, ale také vliv mechanického namáhání. Nejjednodušším způsobem je použít empirické vztahy určující změnu permeability v~závislosti na napěťových (případně deformačních) polích, viz např. \cite{Rutqvist2009Modeling}. Pokročilejší možnost představuje užití víceškálových modelů popsaných v~sekci \ref{sec:hm_homogenizace}. V~rámci námi navržené metodiky použijeme k~popisu změn hydraulické vodivosti kombinaci modelu poroelasticity založeného na Biotově teorii (\cite{Biot1941General}) a Bayesovské inverzní analýzy. 

V této kapitole se věnujeme Biotovu modelu poroelasticity a jeho numerickému řešení. V~sekci \ref{sec:biot} je představen model poroelasticity v~kontinuální oblasti bez puklin. Ten je následně v~sekci \ref{sec:biot_dfm} rozšířen na oblast s diskrétními puklinami a adaptován do kontextu DFM pro popis interakcí v~okolí puklin. Sekce \ref{sec:biot_diskretizace} je věnována diskretizaci a numerickému řešení implementovanému v~SW Flow123d (\cite{Brezina2021Flow123d}).
Při řešení dílčích mechanických problémů s kontaktními podmínkami na puklinách je využíváno algoritmů kvadratického programování, které jsou popsány v~sekci \ref{sec:qp} a implementovány v~knihovně PERMON (\cite{Horak2020PERMON}). Bayesovská inverzní analýza bude představena v~následující kapitole.


\subsection{Biotův model}\label{sec:biot}

Mechanické chování horninového masivu vychází z~rovnice rovnováhy
\eq{\label{eq:el} -\div\vc\sigma + \alpha\varrho g\nabla p = \vc f, }
kde $\vc\sigma$ $[\unit{Pa}]$ je tenzor napětí, $\alpha$ $[-]$ je Biotovo číslo, $\varrho$ $[\unit{kg\cdot m^{-3}}]$ je hustota kapaliny, $g$ $[\unit{m\cdot s^{-2}}]$ je gravitační zrychlení, $p$ $[\unit{m}]$ tlaková výška a $\vc f$ $[\unit{N\cdot m^{-3}}]$ zatížení.
Bilanci hmoty v~masivu vyjadřuje rovnice
\eq{\label{eq:fl} \prtl_t(Sp+\alpha(\div\uu)) + \div\vc v = s, }
kde $S$ $[\unit{m^{-1}}]$ je storativita, $\uu$ $[\unit{m}]$ posunutí, $\vc v$ $[\unit{m\cdot s^{-1}}]$ Darcyovská rychlost proudění a $s$ $[\unit{s^{-1}}]$ objemový zdroj kapaliny.
Jako konstitutivní zákony pro tenzor napětí a rychlost uvažujeme Hookův, resp. Darcyův zákon:
\eq{ \vc\sigma=\tn C\eps(\uu), \qquad \vc v = -\tn K\nabla (p+z). }
Pro izotropní prostředí mají Hookův tenzor a tenzor hydraulické vodivosti tvar:
\eq{ \tn C = 2\mu\tn I_4 + \lambda\tn I\otimes\tn I,\qquad \tn K=k\tn I, }
kde $\mu,\lambda$ $[\unit{Pa}]$ jsou Lamého parametry, $k$ $[\unit{m\cdot s^{-1}}]$ je (skalární) hydraulická vodivost a $\tn I,\tn I_4$ jsou jednotkové tenzory druhého, resp. čtvrtého řádu.
Obecně uvažujeme také případy, kdy vodivost může záviset na mechanických veličinách:
\eq{ \tn K=\tn K(\uu,\vc\sigma), \quad \mbox{resp.}\quad k = k(\uu,\vc\sigma).
\label{eqn:nonlinear_cond}}
Příklady empirických vztahů definujících tyto závislosti zmiňujeme v kap. \ref{sec:tsx_hm} a~\ref{sec:nonlinear_fr_effects}.

Biotův model je doplněn počáteční podmínkou pro tlak
\eq{ p(0,\vc x) = p_0 }
a kombinací Dirichletových a Neumannových okrajových podmínek:
$$ p = p_{in} \mbox{ nebo } \vc v\cdot\vc n = v_{in}, $$
$$ \uu = \uu_{in} \mbox{ nebo } \vc\sigma\vc n = \vc t, $$
se zadanou tlakovou výškou $p_{in}$ $[\unit{m}]$, tokem $v_{in}$ $[\unit{m\cdot s^{-1}}]$, posunutím $\uu_{in}$ $[\unit{m}]$ nebo povrchovým napětím $\vc t$ $[\unit{Pa}]$.

\subsection{Biotův model v~DFM popisu}\label{sec:biot_dfm}
S využitím značení zavedeného v~sekci \ref{sec:vek_pocet_dfm} lze rovnice poroelasticity aplikovat na oblast s diskrétními puklinami.
V oblasti $\Omega_m$ uvažujeme rovnice poroelasticity \eqref{eq:el}-\eqref{eq:fl} s neznámými $p_m$ a $\uu_m$.
Tlak $p_f$ a posunutí $\uu_f$ v~puklinách $\Omega_f$ jsou dány semi-diskretizovanými rovnicemi
\eq{ \label{eq:mixed_dim_problem_el_frac} -\fdiv\,\vc\sigma + \alpha_f\varrho g\fgrad p = \vc f_f, }
\eq{ \label{eq:mixed_dim_problem_fl_frac} \prtl_t\left(S_f p_f + \alpha_f\fdiv\,\uu\right) +\fdiv\,\vc v = s_f, }
kde tenzor napětí a tok v~puklině jsou definovány vztahy
\eq{\label{eq:const_laws_f} \vc\sigma_f = \tn C_f\frac{\fgrad\uu+\fgrad^\top\uu}2,\qquad \vc v_f = -\tn K_f\fgrad p. }
% 
Na puklinách dále uvažujeme spojitost toků, tečných napětí a kontaktní podmínky nepronikání:
\eq{ \vc v_m^\pm\cdot\vnu^\pm = \vc v_f^\pm\cdot\vnu^\pm, }
\eq{ (\vc\sigma_m^\pm\vnu^\pm)_t = (\vc\sigma_f^\pm\vnu^\pm)_t, }
\eq{\label{eq:contact} (\Delta^+\vc\sigma + \Delta^-\vc\sigma)\vnu\cdot\vnu \ge 0, \qquad \delta+(\uu_m^+\cdot\vnu^++\uu_m^-\cdot\vnu^-)\ge 0. }
Zde $\vc v_f^\pm$ a $\vc\sigma_f^\pm$ jsou definovány obdobně jako v~\eqref{eq:const_laws_f}, ale s $\fgrad^\pm$ na místě $\fgrad$.


\subsection{Diskretizace a numerické řešení úloh poroelasticity} \label{sec:biot_diskretizace}
Pro numerické řešení sdruženého HM problému v~SW Flow123d je použita formulace v~posunutí, tlakové výšce a rychlosti. Rovnice elasticity je diskretizována metodou konečných prvků, přičemž pole posunů jsou po částech lineární (P1 prvky). Pro rovnici proudění je použita smíšená hybridní metoda konečných prvků, stejně jako v~sekci \ref{sec:transport_num}.
%která zavádí další proměnnou pro stopu tlaku na stěnách elementů. 
%To umožňuje rozpojení stupňů volnosti pro rychlost a následně efektivní lokální výpočet Schurova doplňku.
Oblast s puklinami je rozdělena na simplexové elementy tak, aby puklinové elementy byly totožné se stěnami přilehlých elementů vyšší dimenze.

Diskretizace sdružené poroelastické úlohy vede na nelineární algebraický problém: Hledáme $(\uu^f,\uu^m)$ tak, aby
\begin{subequations}\label{eq:hm_alg}
\eq{ \tn F(\uu^m)\uu^f = \vc f^f - \tn D\uu^m, }
\eq{ \uu^m = \argmin_{\substack{\vc v^m\\\tn B\vc v^m\le\vc c}}\left(\frac12\tn M\vc v^m-\tn D^\top\uu^f-\vc f^m\right) \cdot \vc v^m. }
\end{subequations}
Zde $\uu^f$ a $\uu^m$ jsou vektory stupňů volnosti pro hydraulickou, resp. mechanickou podúlohu.
V případě, kdy nejsou uvažovány kontaktní podmínky ani závislost hydraulické vodivosti na mechanice, se  problém \eqref{eq:hm_alg} redukuje na soustavu lineárních algebraických rovnic
\eq{\label{eq:hm_alg_lin} \begin{bmatrix}\tn F & \tn D\\-\tn D^\top & \tn M\end{bmatrix}\begin{bmatrix}\uu^f\\\uu^m\end{bmatrix} = \begin{bmatrix}\vc f^f\\\vc f^m\end{bmatrix}. }
Na sdružený problém \eqref{eq:hm_alg} je aplikována metoda fixed-stress splitting, vedoucí k~iteračnímu algoritmu: Nechť je dáno $\uu^f_0$. Pak pro $i=0,1,2,...$ hledáme $(\uu^f_{i+1},\uu^m_i)$ jako řešení následujících úloh:
\eq{\label{eq:mech_alg_split} \min_{\uu^m_i} \left(\frac12\tn M\uu^m_i - \vc f^m - \tn D^\top\uu^f_i\right)\cdot\uu^m_i,\quad \tn B\uu^m_i\le\vc c, }
\eq{\label{eq:flow_alg_split} (\tn F(\uu^m_i) + \beta\tilde{\tn F})\uu^f_{i+1} = \vc f^f + \beta\tilde{\tn F}\uu^f_i - \tn D\uu^m_i, }
kde $\tilde{\tn F}$ je matice hmotnosti pro složku $\uu^\lambda$ vektoru $\uu^f=[(\uu^{\vc v})^\top, (\uu^p)^\top, (\uu^\lambda)^\top]^\top$. V tomto algoritmu jsou podúlohy proudění a mechaniky řešeny odděleně. V~případě lineárního poroelastického problému (bez kontaktních podmínek a nelineární vodivosti) pak lze tento iterační postup reprezentovat pomocí blokově trojúhelníkové soustavy:
\eq{ \begin{bmatrix}\tn F+\beta\tilde{\tn F} & \tn O\\-\tn D^\top & \tn M\end{bmatrix}\begin{bmatrix}\uu^f_i\\\uu^m_i\end{bmatrix} = \begin{bmatrix}\vc f^f\\\vc f^m\end{bmatrix} + \begin{bmatrix}\beta\tilde{\tn F} & -\tn D\\\tn O & \tn O\end{bmatrix}\begin{bmatrix}\uu^f_{i-1}\\\uu^m_{i-1}\end{bmatrix},~i=1,2,... }
Volbou iteračního parametru
$$ \beta\approx \frac{\alpha^2\varrho g}{2(\frac23\mu + \lambda)}, $$
kde $\alpha, \varrho, \mu, \lambda$ jsou výše definované materiálové parametry, je dosaženo poměrně rychlé konvergence k~řešení sdružené poroelastické úlohy.

Algebraická úloha \eqref{eq:flow_alg_split} vzniklá z~hydraulického podproblému je redukována pomocí Schurova doplňku na menší úlohu, která je následně řešena předpodmíněnou metodou sdružených gradientů (viz kap. \ref{sec:transport_num}).
Mechanický podproblém \eqref{eq:mech_alg_split} vede na úlohu kvadratického programování s lineárními omezeními (QP), která je popsána v následující sekci \ref{sec:qp}.

\subsection{QP algoritmy pro numerické řešení kontaktních úloh}
\label{sec:qp}
V této sekci popíšeme numerické řešení kontaktní úlohy \eqref{eq:mech_alg_split}. Úlohu můžeme pro větší přehlednost přeznačit následovně:
\begin{equation}
  \argmin_{\vc u} \frac{1}{2}\vc u^{T}\tn M\vc u-\vc u^{T}\vc f
  \quad\text{s.t.}\quad \tn B \vc u \le \bm{c},
  \label{eq:QPprob}
\end{equation}
kde $\tn M \in \Real^{n\times n}$ je pozitivně definitní matice tuhosti, $\vc f:=\vc f^m + \tn D^\top\uu^f_i \in \Real^{n}$ je vektor zatížení, $\vc u:=\vc u_i^m$ reprezentuje hledaný vektor posunů a lineární nerovnost obsahující $\tn B \in \Real^{m\times n}$ a $\bm{c} \in \Real^{m}$ je diskretizovaná podmínka nepronikání.

V kontaktních úlohách je počet nerovnostních omezení ($m$) obvykle mnohem menší než počet stupňů volnosti ($n$). Splnění nerovnostních omezení u primární úlohy \eqref{eq:QPprob} můžeme standardně vynutit zavedením tzv. Lagrangeových multiplikátorů. Následně pak můžeme definovat duální úlohu
\begin{equation}
  \argmin_{\vc \lambda} \frac{1}{2}\vc \lambda^{T}\tn B\tn M^{-1}\tn B^T\vc \lambda-\vc \lambda^{T}\left(\tn B\tn M^{-1}\vc f - \bm{c}\right)
  \quad\text{s.t.}\quad \vc \lambda \ge \bm l,
  \label{eq:QPdual}
\end{equation}
kde $\bm{\lambda} \in \Real^{m}$ je vektor duálních proměnných reprezentující Lagrangeovy multiplikátory a $\bm{l}$ je známý vektor dolních omezení pro $\bm{\lambda}$, jehož tvar vyplyne z~odvození duální úlohy. Duální úloha má nerovnostní omezení v~jednodušším tvaru, má menší počet neznámých a duální matice $\tn B\tn M^{-1}\tn B^T$ je mnohdy lépe podmíněna než primární matice $\tn M$. Mezi řešeními primární úlohy \eqref{eq:QPprob} a duální úlohy \eqref{eq:QPdual} platí následující vztah 
\begin{displaymath}
    \vc u = \tn M^{-1}\left(\vc f - \tn B^T\vc \lambda\right),
\end{displaymath}
tedy stačí se zabývat pouze řešením duální úlohy a následně dopočítat hledané posuny.

Pro jednodušší popis QP řešiče přeznačíme úlohu \eqref{eq:QPdual} pomocí substitucí $\tn A = \tn B\tn M^{-1}\tn B^T$, $\vc b = \tn B\tn M^{-1}\vc f - \bm{c}$, $\vc x = \vc \lambda$, čímž získáme následující minimalizační problém:
\begin{equation}
  \argmin_{\vc x} \frac{1}{2}\vc x^{T}\tn A\vc x-\vc x^{T}\vc b
  \quad\text{s.t.}\quad \vc x \ge \bm{l}.
  \label{eq:QPprob2}
\end{equation}
Poznamenejme, že matici $\tn A$ typicky nesestavujeme, ale pouze postupně aplikujeme jednotlivé operátory na vektor. Akce inverze $\tn M^{-1}$ je realizována buď pomocí přímého řešiče (obvykle pomocí paralelního řešiče MUMPS, viz \cite{amestoy_fully_2001}) nebo metodou sdružených gradientů.

\begin{algorithm}
	\SetStartEndCondition{ }{}{}%
	\DontPrintSemicolon
	\SetKwProg{Fn}{def}{\string:}{}
	\SetKwFunction{Range}{range}%%
	\SetKw{KwTo}{in}\SetKwFor{For}{for}{\string:}{}%
	\SetKwIF{If}{ElseIf}{Else}{if}{:}{elif}{else:}{}%
	\SetKwFor{While}{while}{:}{fintq}%
	\AlgoDontDisplayBlockMarkers\SetAlgoNoEnd\SetAlgoNoLine%
	%\scriptsize
	\KwIn{$\tn A$, $\vc x^{0} \in \Omega_{feas}$, $\vc b$, $\Gamma > 0$, $\overline{\alpha} \in (0,2||\tn A||^{-1}]$}
	$\vc g = \tn A \vc x^{0} - \vc b$, $\vc p = \vc g^f(\vc x^0)$, $k = 0$\;
	\While{$||\vc g^P||$ není dostatečně malé}{
		\If{$||\vc g^c||^2 \le \Gamma^2 ||\vc g^f||^2$}{
			$\alpha_{f} = \max \{ \alpha: \vc x^k - \alpha \vc p \in \Omega_{feas}\}$\;    
			$\alpha_{cg} = \vc g^T \vc p / \vc p^T \tn A \vc p$\;
			\If{$\alpha_{cg} \le \alpha_{f}$}{
				// krok metody sdružených gradientů\;
				$\vc x^{k+1} = \vc x^{k} -\alpha_{cg}\vc p$\;
				$\vc g = \vc g - \alpha_{cg} \tn A \vc p$\;
				$\beta = \vc p^T \tn A \vc g^f/ \vc p^T \tn A \vc p$\;
				$\vc p = \vc g^f -\beta \vc p$\;
			}\Else{
				// expanzní krok s pevnou délkou\;
				$\vc x^{k+\frac{1}{2}} = \vc x^{k} - \alpha_{f} \vc p$\;
				$\vc g = \vc g - \alpha_{f} \vc p$\;
				$\vc x^{k+1} = P_{\Omega_{feas}} (\vc x^{k+\frac{1}{2}} -\overline{\alpha} \vc g^f)$\;
				$\vc g = \tn A \vc x^{k+1} -\vc b$\;
				$\vc p = \vc g^f$
			}
		}\Else{
			// proporční krok\;
			$\alpha_{cg} = \vc g^T \vc g^c / (\vc g^c)^T \tn A \vc g^c$\;
			$\vc x^{k+1} = \vc x^{k} -\alpha_{cg}\vc g^c$\;
			$\vc g = \vc g - \alpha_{cg} \tn A \vc g^c$\;
			$\vc p = \vc g^f$\;
		}
		$k=k+1$\;
	}
	\KwOut{$\vc x^{k+1}$}
	
	\caption{MPRGP}
	\label{alg:mprgp}
\end{algorithm}

Úlohu \eqref{eq:QPprob2} řešíme pomocí tzv. MPRGP algoritmu (\cite{dostal_optimal_2009}), viz Algoritmus \ref{alg:mprgp}. Tento algoritmus je založen na práci s tzv. aktivní množinou indexů $\cal{A}$ a volnou množinou indexů $\cal{F}$ k vektoru $\bm x$:
\begin{equation}
	\mathcal{A}:=\mathcal{A}(\bm x)=\{j:\; {x}_j = l_j \},\quad
	\mathcal{F}:=\mathcal{F}(\bm x)=\left\{ j:\; {x}_j > l_j \right\}.
\end{equation}
Tyto množiny aktualizujeme při každé změně vektoru $\bm x$, tj. při definování vektorů $\bm x^0$, $\bm x^{k+1}$ nebo $\bm x^{k+1/2}$. Dále v Algoritmu \ref{alg:mprgp} pracujeme s vektorem $\bm{g}$, který reprezentuje gradient kvadratické funkce, jenž minimalizujeme. Jeho komponenty označíme $g_1, g_2,\ldots, g_m$. Při každé aktualizaci gradientu $\bm{g}$ také definujeme tzv. volný (free) gradient $\bm{g}^f$ s komponentami
\begin{flalign*}
  g_j^f = \begin{cases}
  	0 \quad &\text{pokud}\quad j\in\mathcal F,\\
    g_j \quad &\text{jinak},
    \end{cases},\quad j=1,2,\ldots,m,
\end{flalign*}
useknutý (chopped) gradient $\bm{g}^c$ mající komponenty
\begin{flalign*}
  g_j^c =\begin{cases}
  \min(g_j,0) \quad&\text{pokud}\quad  j\in\mathcal A,\\
  0 \quad &\text{jinak},
  \end{cases},\quad j=1,2,\ldots,m,
\end{flalign*}
a projektovaný gradient $\bm{g}^{P}= \bm{g}^f + \bm{g}^c$, jehož norma je přirozeným zastavovacím kritériem algoritmu. Všimněme si, že definice $\bm{g^f}$ a $\bm{g^c}$ jsou založené na práci s~aktivní množinou indexů $\cal{A}$ a volnou množinou indexů $\cal{F}$. Za pomocí volného a useknutého gradientu pracuje Algoritmus \ref{alg:mprgp} zvlášť s komponentami v~aktivní množině $\cal{A}$ a ve volné množině $\cal{F}$. 

V Algoritmu \ref{alg:mprgp} rozlišujeme tři typy kroků: krok metody sdružených gradientů (řádky 8-11), expanzní krok s pevnou délkou (řádky 14-18) a proporční krok (řádky 21-24). Pokud dominuje norma volného gradientu $\bm{g}^f$ vůči normě useknutého gradientu $\bm{g}^c$ (tj. platí podmínka na řádku 3), pak nastává buď krok metody sdružených gradientů nebo expanzní krok s pevnou délkou. První možnost znamená udělat jednu iteraci metody sdružených gradientů za předpokladu, že neopustíme přípustnou množinu $\Omega_{feas}=\{\vc x : \vc x \ge \bm{l}\}$. Pokud by tento krok vedl k~opuštění přípustné množiny, tak místo něj uděláme expanzní krok, který rozšíří aktivní množinu $\cal{A}$. Expanzní krok provádíme ve směru sdružených gradientů z~první možnosti, a to o~maximální délce, při které výsledný vektor zůstane v~přípustné množině. Poté aplikujeme projekci volného gradientu s konstantní délkou kroku $\overline{\alpha} \in (0,2||\tn A||^{-1}]$ (řádek 18), kde projekce na přípustnou množinu $\Omega_{feas}$ je definována jako
\begin{displaymath}
	\left[P_{\Omega_{feas}}(\vc x)\right]_j = \max(l_j,x_j).
\end{displaymath}
Pokud naopak dominuje norma useknutého gradientu $\bm{g}^c$ nad normou volného gradientu $\bm{g}^f$, tak algoritmus udělá proporční krok, který provede posun vektoru $\bm x^k$ v~opačném směru ke $\bm{g}^c$ (viz řádek 22). Díky této akci dojde k~uvolnění indexů z~aktivní množiny $\cal{A}$, která se tak zmenší.

Implementace Algoritmu \ref{alg:mprgp} je k~dispozici v~knihovně PERMON \cite{Horak2020PERMON}, jejíž propojení s knihovnou Flow123d je součástí příslušného softwarového řešení této metodiky. Úlohu \eqref{eq:QPprob} je možné navíc řešit paralelně, kde paralelizace spočívá v~řádkové distribuci matic a vektorů při využití MPI knihovny. Nicméně pro masivní paralelizaci úlohy a zlepšení konvergence, pracujeme na využití popsaných metod spolu s metodou rozložení oblasti typu FETI. Toto rozšíření o~metodu FETI je popsáno v~sekci \ref{sec:feti}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stochastické metody}
\label{sec:nejistoty}

Nedílnou součástí navrhované metodiky je výzkum vlivu nejistot na hodnotu indikátoru bezpečnosti. 
Nejistoty hrají v~geotechnických aplikacích výraznou roli. V~našem případě budeme pro jednoduchost uvažovat nejistoty způsobené nepřesností v~měřeních, v~parametrech uvažovaných matematických modelů nebo v~distribuci puklinové sítě. 

K popisu nejistot budeme používat stochastické metody. V~takovém případě jsou vybrané vstupní parametry náhodnými veličinami. To zjednodušeně řečeno znamená, že nejsou dány konkrétními hodnotami, ale daným rozdělením pravděpodobností pro jednotlivé jejich hodnoty. Označíme-li $X$ jako náhodný vektor vstupních parametrů, potom výstupem, který reprezentuje bezpečnostní indikátor $I(p)$, není konkrétní hodnota, ale také náhodná veličina $Y$ daná svým pravděpodobnostním rozdělením závisejícím na rozdělení $X$, tj. $Y:=Y(X)$. 

Pro známé rozdělení nejistot v parametrech $X$ a deterministickou funkci $Y(X)$ je získání rozdělení $Y$ netriviální úlohou označovanou též jako \emph{propagace nejistot}. V části \ref{sec:mlmc} této kapitoly popíšeme víceúrovňovou metodu Monte Carlo, kterou v kombinaci s metodou maximální entropie používáme pro efektivní aproximaci rozdělení $Y$. 

Zbytek kapitoly je věnován popisu Bayesovské inverzní analýzy, která umožňuje určit pravděpodobnostní rozdělení vstupních parametrů modelu na základě naměřených dat odpovídajících výstupům modelu, ale zatížených nejistotami (šumem), viz sekce \ref{sec:Bayes}. Bayesovskou inverzní analýzu využíváme pro získání parametrů poroelastického modelu vzniku EDZ.

\subsection{Propagace nejistot}
\label{sec:mlmc}

Indikátor bezpečnosti $I(p)$ chápeme jako náhodnou veličinu $Y(X)$ danou deterministickou funkcí $Y$ (transportním modelem) a \sy{pravděpodobnostním} rozdělením jeho parametrů $X$. Jako parametry zatížené nejistotou uvažujeme zejména parametry EDZ získané Bayesovskou inverzí (viz níže), dále pak náhodnou síť puklin a předepsaná rozdělení pro parametry horniny. Pro komplexnější scénáře lze uvažovat též o~zahrnutí náhodného rozdělení poruch jednotlivých ÚOS v čase. Cílem této kapitoly je popsat metodu pro charakterizaci \sy{pravděpodobnostního} rozdělení $Y$, tedy rozdělení indikátoru bezpečnosti.

Z praktických důvodů se omezíme na případ, kdy funkci $Y$ umíme pouze vyhodnotit pro danou realizaci parametrů $X_n$. Základní úlohou je odhad střední hodnoty $\E Y$ veličiny $Y$. Jak později ukážeme, efektivní odhad střední hodnoty je postačující k~aproximaci celého rozdělení. Klasickou metodou odhadu je metoda Monte Carlo, spočívající v jednoduché aproximaci aritmetickým průměrem vzorků:
$$
   \E Y \approx \avg{Y}_N := \frac{1}{N}\sum_{n=1}^N Y(X_n).
$$
Dle centrální limitní věty má rozdíl $\E Y - \avg{Y}_N$ přibližně normální rozdělení s~rozptylem
\begin{equation}
    \label{eq:CLV}
    V(\avg{Y}_N) = \frac{V(Y)}{N}.
\end{equation}
Rozptyl zde můžeme chápat jako kvadrát chyby aproximace střední hodnoty pomocí průměru. Tato chyba tedy klesá s počtem vzorků $N$ velmi pomalu, což znemožňuje získat přesné odhady pro výpočetně náročné funkce $Y$. Podle \eqref{eq:CLV} však zmenšení chyby odhadu můžeme dosáhnout též zmenšením rozptylu $V(Y)$. Toho lze dosáhnout aproximací funkce $Y$ pomocí funkce $\tilde Y$, která bude výpočetně laciná. Odhad pak rozložíme:
$$
    \avg{Y}_{N,M} = \avg{Y - \tilde Y}_N + \avg{\tilde Y}_M = 
    \frac{1}{N}\sum_{n=1}^N Y(X_n) - \tilde Y(X_n) + \frac{1}{M}\sum_{m=1}^M \tilde Y(X_m),
$$
přičemž pro dobrou aproximaci máme malý rozptyl $V(Y - \tilde Y)$, a tudíž pro dosažení dobrého odhadu stačí malé $N$, tj. málo vyhodnocení drahé funkce $Y$. Naopak budeme mít velké množství laciných vyhodnocení $\tilde Y$ dané číslem $M$. 

\subsubsection{Víceúrovňová metoda Monte Carlo}
Zobecněním tohoto postupu pro sekvenci aproximací dospějeme k~odhadu pomocí víceúrovňové metody Monte Carlo (MLMC), \cite{Giles2015Multilevel}.  K detailnímu modelu $Y = Y^L$ sestavíme posloupnost aproximací $Y^l$, $l=1,\dots,L$. Kvůli konzistenci dalších zápisů definujeme také $Y^0 = 0$. Potom pro střední hodnotu veličiny $Y$ platí:
$$
\E Y = \E(Y^L - Y^{L-1}) + \dots + \E(Y^2- Y^1) + \E (Y^1  - Y^0)
$$
Odhadem uvedených středních hodnot pomocí průměrů dostáváme MLMC estimátor:
\begin{equation}
    \avg{Y}_{\vc N} = \sum_{l=1}^L \avg{Y^l - Y^{l-1}}{N_l} =  
    \sum_{l=1}^L \frac{1}{N_l} \sum_{n=1}^{N_l} \big\{Y^l(X_n) - Y^{l-1}(X_n)\big\},
\end{equation}
přičemž $\vc N = (N_L,\dots, N_1)$ je vektor počtu vzorků na jednotlivých úrovních. Dvojice aproximací $(Y^l(X_n), Y^{l-1}(X_n))$ na jedné úrovni musí být vypočtena pro společný vstupní vektor $X_n$. Pro úroveň $l$ označíme rozptyl rozdílu $Y^l - Y^{l-1}$ jako $V_l$ a výpočetní cenu jedné dvojice na úrovni $l$ jako $C_l$. Celkový rozptyl a~cena estimátoru jsou pak určeny následovně:
$$
    V = \sum_{l=1}^L \frac{V_l}{N_l},\qquad C = \sum_{l=1}^L N_l C_l.
$$
Pro daný požadovaný rozptyl $V$ lze najít vektor $\vc N$ minimalizující $C$ ve tvaru $N_l \propto \sqrt{V_l C_l^{-1}}$. Pokud rozptyl $V_l$ dostatečně rychle klesá a $C_l$ dostatečně roste s~rostoucí úrovní $l$, má MLMC estimátor výpočetní náročnost úměrnou náročnosti výpočtu vzorku $Y^L$. V praxi lze dobrých odhadů dosáhnout v časech, které odpovídají zhruba stovce vyhodnocení jemného modelu.  

Pro aplikaci metody MLMC potřebujeme vytvořit k~detailnímu modelu transportu $Y^L$ vhodné hrubší aproximace. Ty zkonstruujeme pomocí numerické homogenizace popsané v kapitole \ref{sec:transport_multiscale}. 


\subsubsection{Aproximace hustoty pravděpodobnosti}
MLMC estimátor poskytuje výpočetně efektivní odhad střední hodnoty. Přirozeně lze pak dostat také odhad rozptylu aplikací estimátoru na odvozenou veličinu 
$(Y - \avg{Y})^2$ a podobně lze pak obdržet i vyšší centrální momenty. Z vypočtených momentů pak lze teoreticky aproximovat hustotu pravděpodobnosti veličiny $Y$. Kvůli špatné podmíněnosti takové úlohy je však lépe odhadnout střední hodnoty zobecněných momentů:
$$
    \mu_r := \avg{\phi_r(Y)}_{\vc N},
$$
kde $\phi_r(x)$, $r=1,\dots, R$, jsou vhodné bázové funkce. Obvykle uvažujeme ortogonální systém polynomů na kompaktním intervalu $\omega$, na kterém budeme hledat aproximaci hustoty. Uvedené odhady $\mu_r$ lze získat pro jednu sadu vzorků a jejich výpočet tak není náročnější než výpočet odhadu jedné střední hodnoty.

Následně hledáme hustotu pravděpodobnosti, která má shodné zobecněné momenty $\mu_r$. Hledáme hustotu s maximální entropií, tedy s minimem další informace.  Řešením příslušné variační úlohy nalezneme hustotu
ve tvaru:
$$
    \tilde \rho(Y) = \exp\Big\{ \sum_{r=0}^{R}\lambda_r \phi_r(Y)\Big\},
$$
přičemž předpokládáme, že $\phi_0(x) = 1$ a $\mu_0 = 1$, což zajišťuje správnou normalizaci.

Vektor $\vc \lambda$ získáme řešením nelineární soustavy rovnic:
$$
    \int_\omega \phi_r(t) \exp\Big\{ \sum_{r=0}^R\lambda_r \phi_r(Y)\Big\} d t = \mu_r.
$$
Uvedený postup je nazýván metoda maximální entropie a popsán v kombinaci s MLMC
v práci \cite{Bierig2016Approximation}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Bayesovská inverzní analýza}
\label{sec:Bayes}

V této sekci popíšeme základní myšlenky a cíle inverzní analýzy, stručně zmíníme deterministický přístup k~inverzní analýze, ale hlavně popíšeme stochastický přístup založený na Bayesovské inverzní analýze a související metody řešení. Aplikaci Bayesovské inverzní analýzy na určení parametrů poroelastického modelu na základě dostupných měřeních v~rámci TSX experimentu popíšeme v~kapitole \ref{sec:TSX}.
 
%V kontextu analýzy EDZ je naším cílem určit co nejlépe vlastnosti horniny v blízkém okolí ražených chodeb na základě dostupných změřených dat. Pro budoucí predikce potřebujeme určit parametry poroelastického modelu, tedy vyřešit inverzní problém. Vzhledem k~omezenému množství dat a velkých nejistot v určení parametrů modelu, jeví se jako vhodný přístup popsat
%tyto parametry jako náhodné veličiny pomocí pravděpodobnostního rozdělení. K tomuto účelu využijeme Bayesovskou inverzní analýzu popsanou níže.

%Ze znalosti pravděpodobnostního popisu parametrů budeme dále těžit v navazujícím transportním modelu, do kterého jako vstupní parametry potřebujeme hydraulickou vodivost a porozitu. Transportní model bude variován prostřednictvím víceúrovňové Monte Carlo metody, kde budeme jednotlivé vzorky vybírat na základě určených pravděpodobnostních rozdělení vybraných parametrů. Pro vytěžení již spočítaných dat vzorky vybírány náhodně z již spočtených poroelastických výpočtů. \sy{[Standa: Prosím o~revizi poslední věty.]}

% TSX měření - dostupná data a další metody

%V této části si představíme úlohu pro nalezení vstupních hodnot parametrů poroelastického
%modelu a základní rozdíl mezi deterministickým a Bayesovským přístupem
%k řešení.

\subsubsection{Obecný popis inverzní úlohy}

Uvažujeme operátor $G:\Real^{N}\rightarrow\Real^{M}$ reprezentující tzv. dopředný model, který vektoru s $N$ vstupními parametry přiřadí vektor s $M$ výstupními hodnotami. V~našem případě bude vstupní vektor obsahovat parametry poroelastického modelu a operátor $G$ bude zahrnovat poroelastický výpočet, přičemž z~jeho řešení sestavíme vhodně zvolený vektor výstupů. Dále předepisujeme vektor naměřených hodnot $\boldsymbol{y}\in\Real^{M}$ získaných na základě experimentálního výzkumu. Inverzní úloha je potom formulovaná následovně: 
\[
\text{najdi }\boldsymbol{u}\in\Real^N:\,\boldsymbol{y}=G\left(\boldsymbol{u}\right),
\]
tedy hledáme takové vstupní parametry, pro které je výstup z~dopředného modelu totožný s naměřenými hodnotami. V~experimentální praxi však matematický model dokonale nekopíruje skutečné chování horninového masivu a navíc každé měření má omezenou přesnost. Tedy realističtější inverzní úloha
bude vypadat následovně: 
\[
\text{najdi }\boldsymbol{u}\in\Real^N:\,\boldsymbol{y}=G\left(\boldsymbol{u}\right)+\boldsymbol{\eta},
\]
kde vektor $\boldsymbol{\eta}$ představuje chybu (modelu a měření). Tedy
naše měření je dáno jako součet výstupu modelu a chyby, přičemž tato chyba
je samozřejmě neznámá. V~tomto případě nemůžeme mluvit o~nalezení přesného řešení inverzní úlohy.

Složitost řešení inverzních úloh obecně závisí na vlastnostech funkce $G$ a~časové náročnosti sestavení tohoto operátoru, která může být značná, pokud je sestavení $G$ závislé na řešení rozsáhlých úloh s parciálními diferenciálními rovnicemi, což je i námi uvažovaný případ.

\subsubsection{Deterministický přístup k inverzní analýze}

Vzhledem k~tomu, že výše uvedené inverzní úlohy nemusí být řešitelné, snažíme se pomocí deterministických metod minimalizovat rozdíl mezi $G\left(\boldsymbol{u}\right)$
a $\boldsymbol{y}$. To lze učinit například pomocí metody nejmenších čtverců, přičemž řešíme následující minimalizační problém
\[
\min_{\boldsymbol{u}\in\Real^{N}}\left\Vert G\left(\boldsymbol{u}\right)-\boldsymbol{y}\right\Vert ^{2}.
\]
Řešením je pak vektor $\boldsymbol{u}$, pro který výraz $\left\Vert G\left(\boldsymbol{u}\right)-\boldsymbol{y}\right\Vert$ nabývá nejmenší hodnoty.
Takové řešení vždy existuje, ale nemusí být unikátní. Navíc, chyba $\boldsymbol{\eta}$ může
způsobit, že některé parametry budou velmi špatně odhadnuty. Je běžné, že
parametry, které ovlivňují výstup z~funkce $G$ jen velmi málo, bývají metodou nejmenších čtverců odhadovány velice nepřesně a vzdáleně od očekávaných hodnot.

Tento problém do jisté míry řeší tzv. regularizace. Regularizace přidává penalizační člen, který se zvětšuje s rostoucí vzdáleností $\boldsymbol{u}$ od apriorně odhadovaného vektoru vstupních hodnot, který označíme $\boldsymbol{u}_{0}$. Vektor $\boldsymbol{u}_{0}$ může být například odhadován z~jiných měření nebo z~expertní zkušenosti. Minimalizace s regularizací má potom tvar
\[
\min_{\boldsymbol{u}\in\Real^{N}}\left\Vert G\left(\boldsymbol{u}\right)-\boldsymbol{y}\right\Vert ^{2}+\rho\left\Vert \boldsymbol{u}-\boldsymbol{u}_{0}\right\Vert ^{2},
\]
kde parametr $\rho>0$ určuje sílu regularizace. Poznamenejme, že může být obtížné vhodně zvolit tento regularizační parametr. 

\subsubsection{Bayesovská inverzní analýza a metody jejího řešení}

Bayesovská inverzní analýza se dá do jisté míry chápat jako rozšíření
deterministického přístupu a bude přinášet dvě hlavní výhody:
\begin{enumerate}
\item není třeba určovat regularizační parametr $\rho$,
\item výsledek není ve formě konkrétního vektoru $\boldsymbol{u}$, nýbrž ve formě pravděpodobnostního
rozdělení obsahujícího všechna možná řešení, včetně jejich pravděpodobnosti.
\end{enumerate}
Naopak, vzhledem k~větší komplexitě výstupu z~Bayesovské inverzní analýzy, musíme počítat s větší časovou náročností oproti deterministickým metodám. 

Pro formulování úlohy Bayesovské inverzní analýzy budeme předpokládat následující:
\begin{enumerate}
\item Uvažujeme aditivní model $\boldsymbol{y}=G\left(\boldsymbol{u}\right)+\boldsymbol{\eta}$,
kde $\boldsymbol{u}$ a $\boldsymbol{\eta}$ jsou nyní náhodné vektory, které jsou na sobě nezávislé.
\item Je dáno apriorní rozdělení $\pi_0$ náhodné veličiny $\boldsymbol{u}$,
tedy veškeré přípustné hodnoty parametrů včetně jejich pravděpodobnosti.
\item Je dáno pravděpodobnostní rozdělení $f_{\boldsymbol{\eta}}$ chyby měření $\boldsymbol{\eta}$.
\end{enumerate}
Apriorní rozdělení $\pi_{0}$ je rozšířením předpokladu očekávané hodnoty $\boldsymbol{u}_{0}$
z deterministického přístupu, $\boldsymbol{u}_{0}$ je obvykle střední
hodnota (případně modus) $\pi_{0}$. Rozdělení chyby $f_{\boldsymbol{\eta}}$
popisuje naši znalost o~přesnosti modelu a měření $\boldsymbol{y}$.
Předpokládáme, že $f_{\boldsymbol{\eta}}$ má nulovou střední hodnotu.

Z pohledu pravděpodobnosti nás nyní zajímá podmíněné pravděpodobnostní
rozdělení $\pi\left(\mathbf{u}|\mathbf{y}\right)$, tedy rozdělení
$\boldsymbol{u}$ při zohlednění měření $\boldsymbol{y}$. Bayesova
věta nám dává následující vztah
\begin{equation}
\begin{array}{c}
\pi\left(\mathbf{u}|\mathbf{y}\right)=\frac{f_{\boldsymbol{\eta}}\left(\mathbf{y}-G\left(\mathbf{u}\right)\right)\pi_{0}\left(\mathbf{u}\right)}{\intop_{\Real^{N}}f_{\boldsymbol{\eta}}\left(\mathbf{y}-G\left(\mathbf{u}\right)\right)\pi_{0}\left(\mathbf{u}\right)\text{d}\mathbf{u}},\end{array}\label{eq:Bayes_normalized-1}
\end{equation}
kde jmenovatel je konstantní a dá se tedy chápat jako normalizační
člen. Tento normalizační člen ovšem nelze v~praktických úlohách přesně vyhodnotit,
a proto budeme pracovat s formou Bayesovy věty
\begin{equation}
\begin{array}{c}
\pi\left(\mathbf{u}|\mathbf{y}\right)\propto \tilde{\pi}\left(\mathbf{u}|\mathbf{y}\right) =  f_{\boldsymbol{\eta}}\left(\mathbf{y}-G\left(\mathbf{u}\right)\right)\pi_{0}\left(\mathbf{u}\right),\end{array}\label{eq:Bayes_theorem-1}
\end{equation}
kde $\propto$ značí úměrnost, tj. platí $\pi\left(\mathbf{u}|\mathbf{y}\right)= q\tilde{\pi}\left(\mathbf{u}|\mathbf{y}\right)$, kde $q>0$ je multiplikativní konstanta. 

Pro popis neznámé hustoty pravděpodobnosti $\pi\left(\mathbf{u}|\mathbf{y}\right)$ budeme využívat vzor\-ky vygenerované pomocí Metropolis-Hastings (MH) algoritmu, viz Algoritmus \ref{alg:MH-2}. Na vstupu předpokládáme hustotu $\tilde\pi\left(\mathbf{u}|\mathbf{y}\right)$ definovanou pomocí \eqref{eq:Bayes_theorem-1} a navíc tzv. návrhovou hustotu $q\left(x|k\right)$. Jedná se o~rozdělení pravděpodobnosti určující návrhy kroků náhodné procházky po parametrickém prostoru. Pro veškeré praktické výpočty budeme používat $q\left(x|k\right)$ jako vícerozměrné normální rozdělení se střední hodnotou rovnou $k$ a diagonální kovarianční maticí. Kovarianční matice je volená v~závislosti na konkrétní úloze, a to s ohledem na apriorní rozdělení, viz \cite{beresova_bayesian_2022}.

\vspace{2ex}
\begin{algorithm}[H]
\textbf{Vstup:} počáteční snímek $k_{1}$, hustota $\tilde{\pi}\left(x|y\right)$ daná vztahem (\ref{eq:Bayes_theorem-1}), návrhová hustota $q\left(x|k\right)$.

Pro $t=1,2,\ldots,T$:
\begin{enumerate}
\item Vygenerujeme snímek $x$ z~rozdělení daného hustotou $q\left(x|k_{t}\right)$.
\item Spočteme pravděpodobnost přijetí
\begin{equation*}
\alpha\left(k_{t},x\right)=\min\left\{ 1,\frac{\tilde\pi\left(x|y\right)q\left(k_{t}|x\right)}{\tilde\pi\left(k_{t}|y\right)q\left(x|k_{t}\right)}\right\}. \label{eq:MH_alpha-2}
\end{equation*}
\item Vygenerujeme prvek $\xi_t$ z~rovnoměrného rozdělení intervalu $[0,1]$. 
\item Pokud $\xi_t>\alpha\left(k_{t},x\right)$ nepřijímáme $x$ a tedy nastavíme $k_{t+1}=k_{t}$.
\item Pokud $\xi_t\leq\alpha\left(k_{t},x\right)$ přijímáme $x$ a tedy nastavíme $k_{t+1}=x$.
\end{enumerate}
\textbf{Výstup:} posloupnost korelovaných vzorků $\left\{ k_{t}\right\} _{t=1}^{T}$ vygenerovaných z~hustoty pravděpodobnosti $\pi\left(\mathbf{u}|\mathbf{y}\right)$

\caption{\label{alg:MH-2}Metropolis-Hastings (MH)}
\end{algorithm}
\vspace{2ex}

Poznamenejme, že Algoritmus \ref{alg:MH-2} je časově velice náročný neboť potřebujeme provést velké množství kroků algoritmu a~v~každém z~nich potřebujeme vyhodnocení funkce $G$. V~našem případě vyhodnocení funkce $G$ představuje řešení parciálních diferenciálních rovnic, konkrétně řešení poroelastické úlohy. 

Pro urychlení výpočtu se dá využít modifikace MH algoritmu, takzvaný MH se zpožděným přijetím (DAMH), viz Algoritmus \ref{alg:DAMH-2}.  Tato modifikace MH algoritmu používá aproximaci neškálované aposteriorní hustoty $\rho \left(x|y\right) \approx \tilde{\pi}\left(x|y\right)$ s~cílem omezit množství vyhodnocení funkce $G$. DAMH vyhodnocuje v~$\tilde{\pi}\left(x|y\right)$ (a~tedy funkci $G$) pouze pro snímky, které označila aproximovaná hustota $\rho \left(x|y\right)$ jako dostatečně kvalitní (tzv. před-přijetí). 

V praktických aplikacích budeme sestavovat $\rho \left(x|y\right)$ za pomocí výpočetně levné aproximace funkce $G$. Tohoto lze například docílit pomocí aproximace polynomem (tzv. stochastická kolokace) nebo pomocí aproximace radiálními bázovými funkcemi. Obě tyto metody aproximace (tzv. surrogate model) jsou detailně popsány v~\cite{beresova_bayesian_2022} a~implementovány ve využitém softwaru surrDAMH, viz \cite{beresova_surrdamh_2022}.

\vspace{2ex}

\begin{algorithm}[H]
\caption{\label{alg:DAMH-2}Delayed acceptance Metropolis-Hastings (DAMH)}
\textbf{Vstup:} počáteční snímek $k_{1}$, hustota $\tilde{\pi}\left(x|y\right)$ daná vztahem (\ref{eq:Bayes_theorem-1}), aproximace hustoty $\rho \left(x|y\right) \approx \tilde{\pi}\left(x|y\right)$, návrhová hustota $q\left(x|k\right)$.

Pro $t=1,2,\ldots,T$:
\begin{enumerate}
\item Vygenerujeme snímek $x$ z~rozdělení daného hustotou $q\left(x|k_{t}\right)$.
\item Spočteme pravděpodobnost před-přijetí
\begin{equation*}
\tilde\alpha\left(k_{t},x\right)=\min\left\{ 1,\frac{\rho\left(x|y\right)q\left(k_{t}|x\right)}{\rho\left(k_{t}|y\right)q\left(x|k_{t}\right)}\right\}. \label{eq:MH_alpha-2b}
\end{equation*}
\item Vygenerujeme prvek $\xi_t$ z~rovnoměrného rozdělení intervalu $[0,1]$. 
\item Pokud $\xi_t>\tilde\alpha\left(k_{t},x\right)$ nepřijímáme $x$ a tedy nastavíme $k_{t+1}=k_{t}$.
\item Pokud $\xi_t\leq\tilde\alpha\left(k_{t},x\right)$ před-přijmeme $x$ a tedy:
    \begin{enumerate}
        \item Spočteme pravděpodobnost přijetí
        \begin{equation*}
        \alpha\left(k_{t},x\right)=\min\left\{ 1,\frac{\rho\left(x|y\right)\tilde\pi\left(k_{t}|y\right)}{\rho\left(k_{t}|y\right)\tilde\pi\left(x|y\right)}\right\}. \label{eq:MH_alpha-2c}
        \end{equation*}
        \item Vygenerujeme prvek $\vartheta_t$ z~rovnoměrného rozdělení intervalu $[0,1]$. 
        \item Pokud $\vartheta_t>\alpha\left(k_{t},x\right)$ nepřijímáme $x$ a tedy nastavíme $k_{t+1}=k_{t}$.
        \item Pokud $\vartheta_t\leq\alpha\left(k_{t},x\right)$ přijímáme $x$ a tedy nastavíme $k_{t+1}=x$.
    \end{enumerate}
\end{enumerate}
\textbf{Výstup:} posloupnost korelovaných vzorků $\left\{ k_{t}\right\} _{t=1}^{T}$ vygenerovaných z~hustoty pravděpodobnosti $\pi\left(\mathbf{u}|\mathbf{y}\right)$

\end{algorithm}
\vspace{2ex}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Software pro realizaci metodiky}
\label{sec:SW}

%\sy{[Standa: Text v této sekci jsem trochu přeuspořádal. Prosím o kontrolu.]}
% JB: OK 

Stochastický výpočet rozdělení indikátoru bezpečnosti $I(p)$ se zahrnutím vlivu EDZ a nejistot v~modelech poroelasticity a transportu je realizován pomocí specializovaného software Endorse, který je samostatným výsledkem projektu. Schéma stochastického výpočtu indikátoru pomocí SW Endorse je uvedeno na \figref{fig:endorse_lib_scheme}. V~této kapitole uvádíme stručný přehled komponent tohoto softwaru, z~nichž některé jsou obecněji použitelnými knihovnami vyvinutými nebo rozšířenými v~rámci řešení projektu. Software Endorse pak zajišťuje konfiguraci celého výpočtu, komunikaci mezi jednotlivými komponentami, přípravu hrubších aproximací transportního modelu pomocí homogenizace parametrů a vlastní výpočet indikátorů z~výsledků transportu.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\textwidth]{\figpath endorse_scheme.pdf}
    \caption{Stochastický výpočet indikátorů bezpečnosti.}
    \label{fig:endorse_lib_scheme}
\end{figure}

Pro Bayesovskou inverzní analýzu a její univerzální použití byla vyvinuta knihovna surrDAMH \cite{beresova_surrdamh_2022}. Klíčovou vlastností knihovny je akcelerace výpočtů pomocí Metropolis-Hastings
algoritmu s odloženou akceptaci a~využitím zástupných (surrogate) modelů, viz Sekce \ref{sec:Bayes}. Knihovnu používáme k~určení parametrů dopředného HM modelu po vzniku EDZ na základě experimentálních dat.

Dopředný HM (poroelastický) model (kapitola \ref{sec:poroelasticity}) a také transportní model (kapitola \ref{sec:transport_edz}) jsou řešeny pomocí simulátoru Flow123d, viz \cite{Brezina2021Flow123d}. Konkrétně jsou tyto problémy řešeny díky konzistentnímu využití DFM popisu, tedy kombinací kontinua a diskrétních puklin. Pomocí HM modelu jsou vypočtena pole vodivosti a porozity v~EDZ. Tato pole spolu s puklinami vstupují do výpočtu transportního modelu.

V~rámci řešení projektu byl vyvinut poroelastický model se zahrnutím kontaktů a nelineárních efektů na puklinách. Kontaktní úlohy vedou na řešení rozsáhlých úloh kvadratického programování, k~čemuž používáme knihovnu PERMON, viz \cite{Horak2020PERMON}, a její řešiče popsané v sekcích \ref{sec:qp} a \ref{sec:feti}.

Vlastní stochastický výpočet indikátorů s propagací nejistot v~parametrech HM modelu a puklinové sítě v~okolí HÚ provádí knihovna MLMC, viz \cite{Spetlik_MLMC_2022}. Tato knihovna realizuje víceúrovňovou metodu Monte Carlo a~rekonstrukci hustot pomocí metody maximální entropie, viz Sekce \ref{sec:mlmc}. 

Dále byla využita knihovna BGEM, viz \cite{brezina_bgem_2022}, pro tvorbu náhodných sítí puklin a navazující tvorbu výpočetních sítí. Knihovny MLMC a BGEM byly vyvinuty v~rámci předchozího projektu a  částečně rozšířeny pro integraci do softwaru Endorse.

Software Endorse a jeho komponenty používají i mnoho dalších otevřených knihoven. Zejména knihovnu PETSC pro paralelní lineární algebru, knihovny MPI a Mpi4py pro paralelní komunikaci nebo síťovací nástroj GMSH a jeho Python rozhraní. Pro vizualizaci je využíván software Paraview, knihovna VTK a její Python rozhraní Pyvista. 

Přenositelnost softwaru Endorse včetně všech závislostí je zajištěna zapouzdřením do Docker resp. Singularity kontejneru. Součástí software je netriviální navázání MPI spojení mezi paralelně běžícími kontejnery a zajištění interakce z~kontejneru s dávkovým plánovačem PBS, který zajišťuje spouštění dílčích úloh.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{Demonstrace metodiky: predikce indikátoru}
\label{sec:transport_numeric_}

V následujících dvou kapitolách demonstrujeme výpočet bezpečnostního indikátoru podle schématu na \figref{fig:endorse_lib_scheme} z~předchozí kapitoly. Tato kapitola je věnována výsledkům
stochastického transportního modelu, následující kapitola pak demonstraci Bayesovské inverze a hydro-mechanického modelu vzniku EDZ.

Pro kompletní aplikaci navržené metodiky bychom pro jednu lokalitu potřebovali znát:
\begin{itemize}
    \item Priory Bayesovské inverze pro poroelastický model, tj. předběžné odhady pro hodnoty a rozptyly hydraulické vodivosti, porozity, počátečního napětí, Youngova modulu a Poissonova čísla. Uvažováno pro neporušenou horninu.
    \item Kontinuální měření pórových tlaků v~okolí ražené chodby.
    \item Heterogenitu EDZ podél úložné chodby.
    \item Regionální pole tlakové výšky (okrajová podmínka proudění).
    \item Data (poloha, směr, velikost, rozevření, transmisivita) významných známých puklin v~blízkosti HÚ.
    \item Statistiky resp. stochastický model menších puklin  v~okolí HÚ.
    \item Difúzní a disperzní parametry v~hornině.
    \item Difúzní parametr ÚOS (kontejner + bentonit).
    \item Předepsaný vývoj koncentrace stopovače v~ÚOS.    
\end{itemize}
\sy{
Jelikož všechna tato data nemáme dostupná na jedné lokalitě, využili jsme pro demonstraci metodiky data z~dostupných zdrojů na více různých lokalitách. I přes snahu použít pokud možno relevatní hodnoty parametrů je nutné výsledky interpretovat v tomto kontextu. Prezentované výsledky tak především dokládají technickou proveditelnost uvažovaných výpočtů a citlivost na vlastnosti EDZ. 
}

\subsection{Geometrie}
\sy{
Pro transportní model je uvažována parametrizovaná geometrie chodeb úložiště pro koncept horizontálního ukládání.
Koncept horizontálního ukládání v kombinaci s mechanickou ražbou byl zvolen zejména kvůli jednoduší geometrii a možnosti zanedbat vliv samotné ražby na EDZ. Použité metody však umožňují i přípravu složitější geometrie pro případ vertikálního ukládání což demonstruje též zahrnutí páteřní chodby a rozrážek do celkové geometrie. }
Geometrie zahrnuje páteřní chodbu a~volitelný lichý počet úložných vrtů s rozrážkami. V~centrálním úložném vrtu jsou voleny polohy ÚOS se zdrojem kontaminace.
Páteřní chodba a rozrážky mají obdélníkový průřez s kruhovým vrchlíkem (viz \figref{fig:tunel_geometry}). Parametry lze pro páteřní chodbu a rozrážky volit nezávisle, avšak byly zvoleny stejné hodnoty: šířka $w=5.5\,\unit{m}$, výška boku $h=4\,\unit{m}$ a poloměr vrchlíku $r=3.5\,\unit{m}$. Počátek souřadné soustavy je zvolen na počvě rozrážky, konkrétně v~bodě napojení úložného vrtu, který pokračuje v~ose rozrážky. Úložný vrt o~poloměru $1.1\,\unit{m}$ se spodní částí dotýká počvy chodby. Osová rozteč úložných vrtů je $20\,\unit{m}$ a délka $300\,\unit{m}$.

\begin{figure}
    \centering    
    \includegraphics[width=\textwidth]{graphics/tunnel_geometry.pdf}
    \caption{Schéma geometrie pro model transportu.}
    \label{fig:tunel_geometry}
\end{figure}
\sy{
Z časových důvodů byly stochastické výpočty prováděny na výřezu z uvedené geometrie tj. na boxu o výšce 60m, šířce 80m a délce 100m ve směru chodby. Software obsahuje funkce pro homogenizaci heterogenit geometrie a transportních parametrů při použití hrubého kroku sítě na části modelu dále od uvažovaného zdroje kontaminace. Aplikace homogenizací je klíčová pro efektivní výpočet na složitějších geometrii avšak v rámci projektu se již nepodařilo tyto techniky aplikovat v rámci komplexního transportního modelu. 
}
\begin{figure}
    \centering    
    \includegraphics[width=0.95\textwidth]{graphics/case_0_better.pdf}
    \includegraphics[width=0.95\textwidth]{graphics/case_4_better.pdf}
    \caption{Dvě uvažovaná uspořádání deteministických puklin. Nahoře varianta A s jednou puklinou protínající přímo ÚOS a jednou navazující puklinou paralelní s úložným vrtem. Dole varianta B s jedinou deterministickou puklinou zasahující přibližně 10m k úložnému vrtu a horizontálně pokrývající většinu oblasti.
    Ve vertikálním řezu osou úložného vrtu jsou zobrzeny zdroje kontemince na staničení 14m, 32m a 62m. Z pole koncentrací v prvních časech simulací je patrné převažující proudění v EDZ zleva doprava.}
    \label{fig:velke_pukliny}
\end{figure}

\subsection{Puklinové sítě}
\sy{
Uvažujeme dvě množiny puklin: množinu velkých puklin, pro kterou předpokládáme přesnou znalost jejich parametrů, a množinu malých puklin, které jsou náhodné. Množina velkých puklin s průměry v~intervalu $(40, 120)$ $\bunit{m}$ byla vygenerována také náhodně s pevným semínkem náhodného generátoru. Tato množina je stejná pro všechny realizace jemné puklinové sítě. Pro níže uvedenou studii 
uvažujeme dvě uspořádíní deterministických puklin. Variantu A se dvěma puklinami z nichž jedna protíná ÚOS a variantu B s jednou puklinou neprotínající úložnou chodbu, viz \figref{fig:velke_pukliny}.  Jednotlivé polohy zdroje kontaminace se liší právě relativní polohou vůči velkým deterministickým puklinám zejména ve vztahu k převažujícímu směru proudění v EDZ. Dále je vygenerována množina menších puklin s průměrem v~intervalu $(5, 40)$ $\bunit{m}$, která se pro jednotlivé realizace stochastického výpočtu liší.
}
Obě množiny používají stejný stochastický model s parametry převzatými ze zprávy SKB 
\cite{Ohman_Site_2010} pro oblast úložiště (HDR -- repository domain) na lokalitě Forsmark. Data byla získána optickým mapováním stěn vrtů a čerpacími testy z~intervalu délky 1~m, který byl posouván s krokem 0.1~m. Tato podrobná data byla použita pro kalibraci transmisivit a velikostí puklin.

Stochastický model se skládá z~Poissonova procesu pro počet puklin v~oblasti, mocninného zákona pro rozdělení velikosti puklin, Fisherova rozdělení pro směr normály a rovnoměrného rozdělení pro polohu pukliny. Pukliny jsou rozděleny do pěti skupin podle převažujících orientací, každá skupina má samostatné parametry.
Data z~lokality Forsmark obsahují též data pro model transmisivit puklin ve tvaru
$T=a r^b$, kde $r$ je rozměr pukliny a $a, b$ jsou empirické parametry. Transmisivita $\bunit{m^2\cdot s^{-1}}$ je hydraulická vodivost přenásobená rozevřením pukliny $T = k \delta$. Použitím kubického zákona bylo z~transmisivit odvozeno následující rozevření puklin:
$$
    \delta = \sqrt[3]{\frac{12 \eta}{\rho g} a r^b}
$$
kde $\eta$ je dynamická viskozita vody, $\rho$ hustota vody a $g$ gravitační zrychlení.



\subsection{Parametry transportního modelu}
\label{sec:transport_param}
V této části shrneme použité hodnoty dalších parametrů transportního modelu.

\paragraph{Parametry advekce.}

Rychlostní pole způsobující advekci je určeno jednak hydraulickou vodivostí a dále předepsaným tlakem, resp. piezometrickou výškou na hranici. 
\sy{
V Tabulce \ref{tab:conductivity} jsou modře vyznačeny použité hodnoty hydraulické vodivosti. Pro srovnání je uvedena řádová vodivost na puklinách. Pro horninu jsou použity hodnoty získané z~experimentu TSX (viz též kapitola \ref{sec:TSX}), zvětšené o 2 řády, tak abychom se řádově přiblížili předpokládaným hodnotám na lokalitách v ČR. I tak je kontrast vodivostí na puklinách a v~hornině je v~tomto případě extrémní. Výpočet v okolí puklin je na hranici strojové přesnosti a bylo nutné upravit přechodové koeficienty mezi puklinou a kontinuem. Při použití konzistentních dat z jedné lokality by byl kontrast pravděpodobně menší. Vodivost na EDZ je zadána na povrchu chodby
a exponenciálně klesá k hodnotě pro neporušenou horninu ve vzdálenosti 1.5m od stěny chodby. Rozsah EDZ je tak uvažován v maximálním rozsahu pozorovaném v nadloží experimentu TSX. Není také proveden přepočet na menší průměr úložné chodby.
V praxi by tedy rozsah EDZ byl patrně menší. 
}

Piezometrická výška na hranici je interpolována z~regionálního hydrologického modelu pro testovací lokalitu Kraví hora. Převažující směr proudění je ve směru úložné chodby, což maximalizuje možné šíření kontaminace podél EDZ.

\begin{table}[!htb]   
    \centering
    %Forsmark -- pukliny $2.21\cdot 10^{-3}$
    % fracture transmisivity - stress relation 
    % https://onlinelibrary.wiley.com/doi/epdf/10.1111/gfl.12089
    \begin{tabular}{lllc}
        \toprule
        místo & $k$ $[\unit{m\cdot s^{-1}}]$ \\ 
        \midrule
        {\color{blue}Forsmark -- pukliny} & {\color{blue} $2.21\cdot 10^{-3}$}\\
        % Surao TZ 334 (Milan Hokr - neverejne, Astra), str. 85
        granit Sm{\aa}land (\" Asp\" o) & $10^{-14}-10^{-13}$ \\
        EDZ v \" Asp\" o (do 1--5 mm )& $10^{-9}-10^{-11}$ \\
        % Surao TZ 303 (Milan Hokr - neverejne, Astra), str. 61
        Kraví hora (prostor HÚ) & $3.6\cdot10^{-10}-5.2\cdot10^{-9}$\\
        EDZ Grimsel (do 2 m) & $2\cdot 10^{-8}$ \\
        % Rutquist 2009
        TSX (neporušená hornina) &  $5.9\cdot 10^{-15}$\\
        TSX (EDZ) &  $10^{-14}-10^{-11}$\\
        {\color{blue} neporušená hornina} & {\color{blue} $3\cdot 10^{-13}$}\\
        {\color{blue} max EDZ} & {\color{blue} $10^{-9}$}\\
        \bottomrule
    \end{tabular}
    \caption{Srovnání hodnot hydraulické vodivosti na různých lokalitách.}
    \label{tab:conductivity}
\end{table}


\paragraph{Difúzní parametry}
Pro bezpečnost HÚ jsou klíčové difúzní procesy, které způsobují ředění počáteční koncentrace a rozptylují kontaminaci do většího objemu. Tento efekt je způsoben jednak chaotickým pohybem molekul (molekulární difúze) a dále nehomogenitou rychlostního pole na menších škálách (hydrodynamická disperze). 

Skutečná míra difúze v~hornině je udávána pomocí efektivního difuzního koeficientu $D_e$ $[\unit{m^2\cdot s^{-1}}]$. Tento koeficient vůči difuznímu koeficientu ve vodě $D_w$ $[\unit{m^2\cdot s^{-1}}]$  dán vztahem:
$$ D_e = \tau\theta D_w, $$
viz \cite{TZ420_2019}, (3.11). Zde tortuozita $\tau$ $[-]$ a porozita $\theta$ $[-]$ 
zohledňují prodloužení difúzní cesty v~pórovém prostoru. Ve zprávě \cite{TZ303_2018}, Tab. 33 lze nalézt následující hodnoty molekulárních difuzivit v~granitové hornině Českého masivu, viz Tab. \ref{tab:radionuclides}.

\begin{table}[!htb]   
    \centering
    \begin{tabular}{c|ccc|c}
    \toprule
        nuklid & \multicolumn{3}{c|}{$D_e~ \bunit{m^2\cdot s^{-1}}$} & $D_w~ \bunit{m^2\cdot s^{-1}}$\\
        & min & střed & max &\\
        \midrule
        I129 & $1.30\cdot10^{-15}$ & $2.60\cdot10^{-15}$ & $2.60\cdot10^{-14}$ & $2.00\cdot10^{-9}$\\
        Cs135 & $2.13\cdot10^{-14}$ & $1.06\cdot10^{-13}$ & $2.13\cdot10^{-13}$ & $2.10\cdot10^{-9}$\\
        Sr90 & $5.01\cdot10^{-14}$ & $1.02\cdot10^{-13}$ & $1.02\cdot10^{-12}$ & $7.91\cdot10^{-10}$\\
        \bottomrule
    \end{tabular}
    \caption{Difuzivita radionuklidů.}
        \label{tab:radionuclides}
\end{table}
Vzhledem k velmi malé vodivosti horniny uvažujeme též minimální difúzi v~hornině, a sice $1.00\cdot10^{-15}$. Pro širší souvislosti a vliv sorpce odkazujeme na zprávu SKB: \cite{ohlsson_diffusion_1997}.

Hydrodynamická disperze je způsobená heterogenitami hydraulické vodivosti na menších škálách,
a je tedy stejně jako hydraulická vodivost závislá na škále problému.
Monografie \cite{bear_flow_2012} uvádí pro oblasti okolo 10 m podélnou disperzi cca 1 m. 
Pro příčnou disperzi se obvykle uvažuje hodnota řádově nižší, tedy 0.1 m. 
    
Porozita udává objem pórového prostoru (v krystalických horninách realizovaného sítí puklin)
v jednotkovém objemu. Hodnoty porozity se pohybují v~řádu $10^{-4}$--$10^{-3}$.
Konkrétně, v~horninách Českého masivu se pohybuje v~intervalu 0.0014--0.0129, doporučený medián je $0.004$ (viz \cite{TZ303_2018}, kap. B.3.1.2). Podobně, \cite{Rutqvist2009Modeling} uvádí pro URL v~Kanadě hodnotu 0.007. Uvažujeme hodnotu $\theta=0.004$ pro horninu a hodnotu $\theta=1$ na puklinách. Změny pórového prostoru v~EDZ, které způsobují zvýšení vodivosti zvyšují i porozitu, avšak tato změna je velmi malá. 


Únik kontaminace ze zdrojového ÚOS je popsán Robinovou okrajovou podmínkou \eqref{eq:source_bc}. Difúzní parametr $\alpha_B$ reprezentuje difúzi skrze 1~m bentonitové obálky. Uvažujeme zdánlivý difúzní koeficient $D_a=2\cdot 10^{-10}$ zohledňující též retardaci sorpcí. Dle zprávy \cite{Vokal_2010}, Tab. C.2 -- 29 je zvolena maximální hodnota pro nesorbující jód v bentonitu {\it Rokle}. Pro vývoj koncentrace uvnitř UOS se nám nepodařilo dohledat vhodný zdroj. Uvažujeme exponenciální pokles koncentrace z~$0.1\,\unit{g\cdot m^{-3}}$ v~čase 0 o~tři řády na koncentraci $10^{-4}\,\unit{g\cdot m^{-3}}$ v čase 1000 let. Je uvažován čas od porušení UOS.



    
    %     \begin{tabular}{|lllc|}
%         \hline
%         parameter & symbol & value & unit \\
%         \hline
%         hydraulic conductivity & $K$ & $5.886\cdot10^{-15}$ & $\unit{m\cdot s^{-1}}$ \\
%         storativity & $S$ & $7.291\cdot10^{-8}$ & $\unit{m^{-1}}$ \\
% %                 \hfill $S=\rho g(\beta_s + \vartheta \beta_w)$ \\
% %           \begin{itemize}
% %             \item kompresibilita horniny $\beta_s = 0$ $\mathrm{[Pa^{-1}]}$
% %             \item kompresibilita vody $\beta_w = 4\cdot 10^{-10}$ $\mathrm{[Pa^{-1}]}$
% %             \item porozita $\vartheta = 0.007$ $\mathrm{[-]}$
% %           \end{itemize}
%     %
%         initial pressure head & $h_0$ & $300$ & $\unit{m}$ \\
%         Biot-Willis coeffiecient & $\alpha$ & $0.2$ & $\unit{-}$ \\
%         Young modulus & $E$ & $6\cdot10^{10}$ & $\unit{Pa}$ \\
%         Poisson ratio & $\nu$ & $0.2$ & $\unit{-}$ \\
%     \hline
%     \end{tabular}
%\end{table}

% Update on reference model report
% Super container
% - diameter 2.1m
% - length 6.5m
% - bentonite thickness 0.7m
% - different number of 0.5m distance blocks between SC (1, 9, 17 ?? reason)
% some values for diffusioncoeffs, effective porosity, tortuosity and pertition % coefs for essential radionuclides in bentonite
% can not be directly used for granite, with possibly different tortuosity
%
% References:
% (diffusion coefs.)
% [83]
% [18]
% [28]
% [29]
%
% partition numbers:
% [30]
% [31]
% [32]

% equivalen flow concept 
% [99]
% [126]  Nykyri M., Norman H., Narcis N., Löfman J., Poteri A., Hautojärvi A.,
%    Radionuclide Release and Transport RTNT-2008, Posiva 2008-06, December 2008



% Bentonite

%  \jb{TODO: přehled řádových hodnot klíčových parametrů modelů pro granit na lokalitách: TSX, Bukov, zprávy SKB Forsmark, další lokality...? Cíl: využit to jako podklad k~parametrům dílčích numerických studiií, získat širší obraz kde se s volbou parametrů pohybujeme}

% Pro žuly na různých lokalitách:
% \begin{itemize}
%     \item molekulární difúze hlavních kontaminantů: (I, Cs, Sr)
%     \item disperze (podélná, příčná)
%     \item porozita
%     \item Youngův modul
%     \item počáteční napětí
% \end{itemize}



%  Speciální kapitolou je předpis zdrojového členu z bentonitu. Podle reportu od Petra Rálka, máme (v CSV) průběh toku na rozhraní bentonit/žula v $[kg/rok]$. Zde je problém, že je časová osa nekonzistentně přeškálovaná na sekundy (asi).  
 
%  Pro reálné izotopy však nestudujeme jejich molekulu a tím příslušnou molární hmotnost, asi je pak lepší počítat s koncentrací v $[mol/m3]$ nebo uvažovat nějakou modelovou molární hmotnost. Jirka  Landa v tabulce nějaké molární hmotnosti pro izotopy má, snad to tedy jsou pro jějakou typickou geochemii. Je to v ARP to by mělo jít dohledat abylo by to vhodné citovat.

%  Pro konzistenci bychom měli počítat s koncentrací Bq/m3. Asi by šlo něco vyčíst z nových zpráv POSIVY.
 
%  \jb{Podle reportu: Working Report 2013-29,
% A Review of Porosity and Diffusion
% in Bentonite, C.A.J. Appelo, 2013; je apparent diffusion (započtena retardase kvůli sorpci a pod) pro Sr v bentonitu $1e-11$ $[m^2/s]$ pro formulaci s koncentrací v jednotce $[mol/m^3]$ ale pro Na je to $1e-10$. }

\subsection{Výpočet bezpečnostních indikátorů}
V této kapitole diskutujeme řešení stochastického transportu pro varianty A a B velkých puklin. Pro tři varianty zdroje kontaminace z ÚOS na staničení 14m, 32m, a 62m. A pro dvě varianty hydraulické vodivosti na EDZ: varianta s EDZ a exponenciálně klesající vodivostí se vzdáleností od stěny úložné chodby (viz kapitola \ref{sec:transport_param}) a varianta bez EDZ s vodivostí neporušené horniny.  Celkem tedy bylo uvažováno 12 variant stochastického výpočtu.

Jedna realizace pro náhodnou síť malých puklin zahrnuje přípravu výpočetní sítě, výpočet prostorově variabilních hodnot porozity a vodivosti, výpočet proudění, výpočet transportu, stanovení hodnoty indikátoru bezpečnosti. 
Pro výpočet je použita síť s krokem 0.3m na stěně úložného vrtu, krokem 1m na vnější hranici EDZ a maximálním krokem 10m na puklinách. Na hranici oblasti je povolen krok až 50m. Výpočetní síť má přibližně 150 tisíc elementů v závislosti na konkrétním uspořádání malých puklin (zhruba 40 puklin). Puklinová síť není přímo regularizována, což vede k tvorbě různých drobných elementů ve výpočetní síti generované knihovnou GMSH. Je proto dále aplikován algoritmus, který regularizuje přímo výpočetní síť i za cenu mírné deformace puklin.   

Výpočet proudění a transportu byl proveden pomocí simulátoru Flow123d. Vývoj pole koncentrace pro jednu realizaci je na Obr. \ref{fig:transport_frames}. Z počátku je patrná významná advekce podél EDZ, zhruba od času 2000 let se kontaminace objevuje na přilehlých puklinách. V čase 10 000 let již je významný transport kontaminace puklinami až k hranici. V koncovém čase je prostorové maximum koncentrace na okraji oblasti což odpovídá sestupující části průnikové křivky. 

\begin{figure}
    \centering    
    % \includegraphics[width=\textwidth]{graphics/transport_frames.pdf}
    \includegraphics[width=\textwidth]{graphics/transport/edz002-s0-t120-fine.png}
    \includegraphics[width=\textwidth]{graphics/transport/edz002-s0-t2300-fine.png}
    \includegraphics[width=\textwidth]{graphics/transport/edz002-s0-t10000-fine.png}
    \includegraphics[width=\textwidth]{graphics/transport/edz002-s0-t200000-fine.png}
    \caption{Vývoj koncentrace kontaminace ve svislém a vodorovném řezu tunelem.}
    \label{fig:transport_frames}
\end{figure}

Po výpočtu simulace transportu je diskrétní pole koncentrace ve výstupních časech vyhodnoceno na 
horizontálních rovinách v blízkosti horní a spodní hranice výpočetní oblasti (viz definice indikátoru v části \ref{sec:transport_indikator}). Z vypočtených hodnot je pro každý výstupní čas spočten kvantil $1-q$ na hladině \sy{významnosti} $q$.
% pro výstupní časy proveden řez polem koncentrace
Na Obr. \ref{fig:indicators_time} jsou typické časové průběhy prostorových kvantilů (tedy funkcí $c_p(t)$, viz sekce \ref{sec:transport_indikator}) pro několik hladin významnosti. V levém sloupci jsou relativně typické průnikové křivky. Horní graf odpovídá případu se špatným propojením puklin, takže šíření k rovinám definice indikátoru trvá déle nežli koncový čas výpočtu. Tento typ průběhu je pochopitelně častý ve variantách výpočtu bez EDZ, nicméně zobrazený průběh je realitvně ojedinělý případ tohoto průběhu pro variantu s EDZ. To lze pozorovat též na Obr. \ref{fig:samples_indicator_vs_time}, kde jsou vyneseny hodnoty indikátorů (časových maxim) spolu s časy jejich dosažení.   

Na pravé straně Obr. \ref{fig:indicators_time} jsou průběhy kvantilů pro případ bez EDZ a variantu~A, kde puklina přímo protíná úložnou pozici na staničení 32m. Při dobrém napojení velké pukliny malými puklinami se kontaminace dostane velmi rychle na definiční rovinu indikátoru, následně dojde k poklesu zdrojového členu čímž dojde i k poklesu koncentrace na této rychlé preferenční cestě. Tato preferenční cesta však díky disperzi stále obsahuje velké množství látky, která je pomalejším příčným prouděním na velké puklině přesunuta k jiné preferenční cestě, kde následně nastane druhé maximum časového průběhu. Pro variantu s EDZ není tento efekt tak výrazný, jelikož látka je zároveň významě transportována skrze EDZ, tedy pryč od pukliny protínající zdroj na staničení 32m. Naopak lze na Obr. \ref{fig:samples_indicator_vs_time}  pozorovat podobný jev pro zdroj na staničení 14m, který je proti směru proudění od pukliny.  


Z kontroly průběhů indikátrů pro větší množství vzorků lze tvrdit, že se podařilo značně stabilizovat výpočet transportu a maximum se výrazně neodchyluje od kvantilů na striktně kladné hladině významnosti. V dalších výsledcích používáme kvantil $1 - 5\cdot 10^{-4}$.


\begin{figure}
    \centering    
    \includegraphics[width=\textwidth]{graphics/indicators_selected.pdf}
    \caption{Ukázky časového vývoje vybraných kvantilů z~hodnot koncentrací na definičních rovinách indikátrou bezpečnoati (viz sekce \ref{sec:transport_indikator}. Hodnota indikátoru (maximum přes časový interval) vyznačena hvězdičkou.}
    \label{fig:indicators_time}
\end{figure}


\begin{figure}
    \centering    
    % \includegraphics[width=\textwidth]{graphics/transport_frames.pdf}
    \includegraphics[width=0.48\textwidth]{graphics/time_conc_A.pdf}
    \includegraphics[width=0.48\textwidth]{graphics/time_conc_B.pdf}
    \caption{Hodnoty indikátorů (desítkový logaritmus) pro nasbírané vzorky spolu s časem jejich dosažení. Vývoj koncentrace kontaminace ve svislém a vodorovném řezu tunelem.}
    \label{fig:samples_indicator_vs_time}
\end{figure}


Z časových důvodů se nepodařilo začlenit algoritmy numerické homogenizace do celkového výpočtu transportu což také znemožňuje demonstraci víceúrovňové metody Monte Carlo. Prezentované výsledky byly tedy vypočteny klasickou metodou Monte Carlo. Výpočet 100 vzorků ve 12-ti variantách trval celkem $9$ hodin při využití $40$ výpočetních jader, tedy přibližně $45$ minut na jeden vzorek.

Obrázek \ref{fig:samples_indicator_vs_time} reprezentuje hodnoty indikátorů a časy jejich dosažení pro sledované varianty. Mimo dříve zmíněných pozorování, lze též vidět zejména vyšší a dřívější hodnoty indikátorů pro případ s EDZ. Značná část bodů je koncentrována podél horní linie, která patně představuje realizaci preferenční cesty pomocí malých puklin. Bez EDZ se tyto preferenční cesty uplatní daleko méně. 
Ve varianě B, kde se velké puklina pouze přibližuje k úložné chodbě je vliv EDZ dokonce několik řádů. 

Častá realizace preferenční cesty pomocí puklin může být způsobena způsobem sestavení stochastického modelu puklin na lokalitě Forsmark. 
Tento stochastický model a závislost transmisivity puklin na jejich velikosti byl kalibrován pomocí DFN modelu proudění. DFN modely z principu vyžadují propojenost puklinové sítě a mohou nadhodnocovat transmisivitu puklin, aby kompenzovali toky skrze neporušenou horninu.

Na Obrázku \ref{fig:samples_indicator_vs_time} dále pozorujeme také poměrně pomalý pokles hodnoty indikátoru (koncentrace) s rostoucím časem jeho dosažení, zde však je pokles způsoben pouze ředěním pomocí difúze a disperze, při zahrnutí sorpcí a rozpadů by byl efekt patrně významější. 


\begin{figure}
    \centering    
    % \includegraphics[width=\textwidth]{graphics/transport_frames.pdf}
    \includegraphics[width=0.9\textwidth]{graphics/mc_cases_A.pdf}
    \includegraphics[width=0.9\textwidth]{graphics/mc_cases_B.pdf}
    \caption{Porovnání aproximace hustot pravděpodobnosti pro hodnoty indikátoru bezpečnosti pro uvažované verianty výpočtu. Vodorovné čárkované čáry vyznačují výběrový medián a kvartily.}
    \label{fig:variant_comparison}
\end{figure}

Na Obrázku \ref{fig:variant_comparison} je uvedena aproximace hustoty pravděpodobnosti indikátoru bezpečnosti pro všechny uvažované varianty stochastického výpočtu. Kvůli použití klasické metody Monte Carlo je použita jednodušší metoda aproximace pomocí jádrových funkcí umístěných v nasbíraných vzorcích. 


Pro variantu A i pro variantu B a pro všechny polohy zdroje pozorujeme výrazný rozdíl mezi případem s EDZ a případem bez EDZ, tedy významný vliv hydraulické vodivosti EDZ na hodnotu indikátoru. Z hlediska bezpečnosti nás bude zajímat zejména porovnání horních kvanitlů získaných rozdělení. Pro horní kvartil pozorujeme rozdíl 1 -- 2 řády pro variantu A puklinou protínající ÚOS. Pro variantu B s velkou puklinou ve větší vzdálenosti je rozdíl 2 -- 4 řády.  Pro případ s EDZ jsou hodnoty indikátoru poměrně vyrovnané pro různé zdroje kontaminace. Rozdíly mezi polohami zdroje jsou však patrné pro variantu bez EDZ. 

Na základě výše uvedených výsledků lze vyslovit hypotézu, že hodnota hydraulické vodivosti EDZ funguje jako filtr vlivu malých puklin. Tedy pro malou hydraulickou vodivost EDZ má na hodnotu indikátoru vliv pouze přítomnost velkých puklin a s klesající hydraulickou vodivostí EDZ roste vliv puklin stále menších velikostí. S rostoucí hydraulickou vodivostí EDZ tak roste počet puklin, které se mohou podílet na tvorbě preferenční cesty a tím roste pravděpodobnost její realizace.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{\sy{TSX experiment a jeho numerická simulace pomocí poroelasticity a Bayesovské inverzní analýzy}}
\section{TSX experiment a Bayesovská inverzní analýza}
\label{sec:TSX}

Transportní parametry EDZ (hydraulické vodivosti, porozity, difúze, disperze a sorpčních koeficienty) mají zásadní vliv na bezpečnost jednotlivých úložných pozic, kterou můžeme kvantifikovat pomocí navrženého indikátoru $I(p)$. Pro tyto parametry chceme proto co nejpřesněji určit
jejich hodnoty včetně popisu jejich heterogenity, zejména v podélném směru úložného vrtu. Pro stochastický výpočet však potřebujeme také popsat nejistoty ve formě rozdělení uvedených parametrů. 

Uvedené transportních parametry jsou závislé na fraktálním charakteru puklinové sítě v hornině, což způsobuje jejich hodnot na měřítku úlohy. Je proto vhodné tyto parametry určovat metodami, které mají integrační objem srovnatelný s měřítkem příslušného transportního modelu. Z podstaty se pak jedná o~inverzní metody. Konkrétně chceme využít Bayesovskou inverzní analýzu popsanou v kapitole \ref{sec:Bayes} pro získání rozdělení uvedených parametrů. 

Pro charakteristiku heterogenity EDZ v rozsahu reálného úložiště by bylo vhodné použití relativně laciných neinvazních geofyzikálních metod. Ovšem tyto metody nejsou schopny poskytnout absolutní hodnoty transportních parametrů, které je nutné získat nezávisle. Pro kalibraci transportních parametrů, zejména hydraulické vodivosti a porozity, jsme proto navrhli využít kontinuální měření pórových tlaků v blízkosti raženého tunelu v kombinaci s poroelastickým modelem vzniku EDZ a Bayesovskou inverzní analýzou. 

V~této kapitole popíšeme určení hydraulické vodivosti a porozity v EDZ z~kontinuálního měření pórových tlaků během ražby a to konkrétně s využitím dat z~experimentu TSX. Nejprve popíšeme vlastní experiment, dále představíme uvažovaný poroelastický model a prezentujeme naše výsledky získané pomocí Bayesovské inverzní analýzy. 

\subsection{Popis TSX experimentu}
\label{subsec:TXS}

Kanadský experiment TSX (Tunnel Sealing Experiment), viz \cite{Chandler2002five}, probíhal v~letech 1995--2002 a byl zaměřen především na výzkum utěsnění chodeb pro hlubinné úložiště, v~rámci kterého probíhalo kontinuální měření pórových tlaků. My se budeme věnovat jedné z~prvních etap tohoto experimentu (po dobu jednoho roku v~rozmezí 1996--1997), kdy se tunel vyrazil. Tato etapa byla unikátní v~tom, že už během ražby probíhala řada měření, na základě kterých se charakterizovaly vlastností horniny (granitu) a EDZ, a byl pozorován vliv EDZ na hydraulické vlastnosti dotčeného masivu.

\begin{figure}[htb!]
    \centering
    \begin{subfigure}[t]{0.35\textwidth}
      \includegraphics[width=\textwidth]{\figpath TSX_situation.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \begin{subfigure}[t]{0.6\textwidth}
      \includegraphics[width=\textwidth]{\figpath TSX_situation2.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \caption{Situační schéma TSX experimentu. Vlevo - schéma podzemních chodeb v~místě podzemní laboratoře. Vpravo - bližší pohled na umístění raženého tunelu (červeně) a okolních chodeb, ze kterých se provádělo měření (zeleně).}
    \label{fig:tsx_situation_scheme}
\end{figure}

Prostorové uspořádání experimentu provedeného v~hloubce okolo 500 metrů přibližuje \figref{fig:tsx_situation_scheme}. Před samotnou ražbou červeně vyznačeného tunelu (o~délce 12~m a průměru 4.375~m) byly na boku a v~nadloží vyraženy pomocné paralelní chodby (zeleně vyznačené) ve vzdálenosti 12~m od budoucího tunelu. Z~boční chodby byl veden horizontální vrt a~z~nadložní chodby svislý vrt směrem k~budoucímu tunelu. Tyto vrty, kolmé k~ose budoucího tunelu, končily 1m od jeho plánované stěny. Byly do nich umístěny 4 senzory pro kontinuální měření hydrostatického tlaku (HGT1-5, HGT1-4, HGT2-4, HGT2-3). Pozice těchto senzorů jsou znázorněny ve 2D řezu na \figref{fig:TSX} nebo \figref{fig:model_01_geometry} -- vertikální vrt (HGT1-5, HGT1-4, označme je pro naše účely V1 a V2) a horizontální (HGT2-4, HGT2-3, označme je pro naše účely H1 a H2).

Dále byly před ražbou změřeny počáteční hodnoty pórového tlaku vody (3~MPa) a tenzoru napětí, přičemž jeho hlavní složky jsou: 60~MPa ve směru testovací komory, 45~MPa v~horizontálním směru kolmém k~ose tunelu a 11~MPa ve vertikálním směru.
Ve skutečnosti byly hlavní směry napětí odkloněny od výše uvedených os o~$8^\circ$, což budeme v~následující numerické simulaci zanedbávat, podobně jako v~\cite{Rutqvist2009Modeling}.

Následně proběhla ražba testovací komory metodou \uv{smooth drill and blast}. Snímání probíhalo kontinuálně během ražby testovacího tunelu i po ní. Naměřený průběh tlaků v~pozorovacích bodech V1, V2, H1 a H2 po dobu jednoho roku je uveden na \figref{fig:TSX}. Vlivem ražby došlo k~výraznému nárůstu tlaku nad tunelem a naopak k~poklesu tlaku na boku tunelu konzistentně s~předpokládanými změnami napětí v~hornině v~důsledku ražby. Relativně velká změna pórových tlaků v~důsledku ražby je způsobena velkými počátečními napětími na lokalitě. Vzhledem k~velmi malé hydraulické vodivosti neporušené horniny (až $10^{-15}\,\unit{m\cdot s^{-1}}$) docházelo následně k~relaxaci tlaků velmi pomalu, zvýšený tlak v~nadloží byl pozorován i~6~let po ražbě. Výjimku ale představuje bod H2, kde dochází k~nečekaně rychlému ustálení na počáteční tlak, který nekoresponduje s~nízkými hodnotami hydraulické vodivosti.

\begin{figure}[!ht]
    \centering
    \includegraphics[width=0.75\textwidth]{graphics/Rutquist_graph_model.pdf}
    \caption{Průběh měřeného a simulovaného pórového tlaku v~okolí nově vyražené chodby, zdroj \cite{Rutqvist2009Modeling}.
    Pozorovací body v nadloží V1, V2 (v originále HGT1-5, HGT1-4) a na boku H1, H2 (HGT2-4, HGT2-3).}
    \label{fig:TSX}
\end{figure}

\cite{Rutqvist2009Modeling} následně provedli numerickou validaci těchto experimentálních výsledků pomocí poroelastických modelů. Pro dobré zachycení průběhu tlaků v~nadloží (V1, V2) museli ale nakonec použít dvě různé sady parametrů, které jsou v \figref{fig:TSX} označeny \uv{Simulation a} a \uv{Simulation b}. Aproximace tlaků v~horizontálním vrtu jim vychází výrazně hůře -- v~bodě H1 se v~numerickém modelu objevují záporné tlaky; ve vzdálenějším bodě H2 numerické modely nedokáží zachytit nečekaně rychlé ustálení počáteční tlaků.
 
Další měření v rámci TSX experimentu, která umožnila charakterizovat hydraulickou vodivost v~okolí testovací komory, zahrnovala  profilování ve vrtech metodami MVP a SEPPI. První metoda, \uv{microvelocity probe}, byla realizována pomocí kompaktních sond s jedním akustickým zdrojem a dvěma geofony. Pomocí sondy bylo provedeno profilování ve vrtech vedených na různé strany od  stěny testovací komory. Metoda SEPPI (viz \cite{bauer_situ_1995}) používá sondu dlouhou 226~cm, zahrnující injektážní komoru o~délce 5~cm, dvojici mikropakrů, dvojici sběrných komor a~dvojici konvenčních pakrů. Malá velikost injektážní komory a velikost mikropakrů umožňuje profilování hydraulické vodivosti s~poměrně detailním prostorovým rozlišením.  Kvalitativně byla pozorována velmi dobrá korelace mezi naměřenou vlnovou rychlostí (MVP) a~hydraulickou vodivostí (SEPPI). Hydraulická vodivost klesala z~hodnot na stěně díla $10^{-12}\,\unit{ m\cdot s^{-1}}$ (bok tunelu), respektive $10^{-10}\,\unit{ m\cdot s^{-1}}$ (nadloží), k~hodnotě $10^{-13}\,\unit{ m\cdot s^{-1}}$ neporušené horniny ve vzdálenosti 1~m od stěny.
Pozdější krátkodobý laboratorní experiment \cite{Souley2001Damageinduced} (opět SEPPI metodou) podpořený i~numerickým modelem poškození ukazuje na hodnoty vodivosti měřeného vzorku o~řád nižší. \cite{Rutqvist2009Modeling} pak uvádí, že kalibrací jejich numerického modelu dospěli k vodivostem neporušené horniny v okolí tunelu menším než $10^{-14}\,\unit{ m\cdot s^{-1}}$.

% TODO: vysvetlit (PE: zatim nechavam zakomentatovano, nez to interne doladime)
% TODO: umime to porovnat v ramci dopredneho modelu s zvysenim vodivosti blizko tunelu? (jak vychazi best-fit??)
%\pe{PE: Standa zde cituje hodnoty vodivosti zrejme z bauer\_situ\_1995 (vodivost je cca $10^7\times$ permeabilita). Jenze Rutqvist2009Modeling popisuje, ze permeabilita z jinych mereni  byla o~1--2 rady vyse a v~modelech ji snizovali radove (2-4 rady), aby  meli rozumny fit -- i proto nase prevzata vodivost $10^{-15}$.}
%\sy{[Standa: Já jsem pouze vycházel z Vašich podkladů, které jste doplnili. Žádné vstupy jsem sám od sebe nedefinoval.]}
% PE: Tak jsem to urcite nemyslel - jen je potreba zde toto diskutovat a vysvetlit, proc potom uvazujeme vodivost o 3 rady jinde.
% Mozna to bude chtit procist vic ty zdroje z Rutquista, coz jsem ted nestihal.

V práci \cite{Rutqvist2009Modeling} jsou experimentální data z kontinuálního měření i data z měření SEPPI použita pro porovnání různých poroelastických modelů se zahrnutím změny permeability. Je při tom použito fitování některých materiálových parametrů na základě měření. V~následujících sekcích na tento výzkum navážeme a pokusíme se o~charakterizaci materiálových parametrů ve smyslu Bayesovské inverzní analýzy.

\subsection{Popis poroelastického modelu}
\label{sec:tsx_hm}


Vycházíme z kontinuálního poroelastického modelu popsaného v~kapitole \ref{sec:poroelasticity} a~implementovaného ve Flow123d. Výpočetní doménou je 2D řez okolím tunelu, viz \figref{fig:model_01_geometry}. Na vnější hranici uvažujeme Dirichletovu okrajovou podmínku pro tlak vody i~posunutí. Ražbu tunelu simulujeme pomocí okrajové podmínky na stěně tunelu, kdy se během prvních 17 dní lineárně snižuje tlak vody k~nule a zvyšuje se normálová síla působící na stěnu tunelu až do hodnoty odpovídající počátečnímu napětí.
Počáteční podmínka zahrnuje pouze tlakovou výšku $300\,\unit{m}$, tj. tlak $3\,\unit{MPa}$.

%\sy{[Standa: Zapoznámkoval jsem níže diskutované porovnání s Kraví horou a či PVP Bukov, které mi přijde  nevhodné a matoucí v této části textu. Doporučil bych to vymazat či %přesunout jinam.]}

% \begin{table}[!htb]
%     \centering
%     \begin{tabular}{lllc}
%         % TODO: sloucit s Tab. 3
%         % sloupecky "stredni hodnota", "odchylka", "rozdeleni" - popsat v textu transformaci na rozdeleni
%         % logaritmus hodnoty je stredni hodnota lognormalniho rozdeleni
%         % sigma v dekadickem logaritmu, samostatny sloupecek
%         % biot: sigma bude odchylka ekvivalentniho normalniho rozdeleni
%         \toprule
%         parametr & symbol & hodnota & jednotka \\
%         \midrule
%         hydraulická vodivost & $k$ & $5.89\cdot10^{-15}$ & $\unit{m\cdot s^{-1}}$ \\
%         storativita & $S$ & $7.29\cdot10^{-8}$ & $\unit{m^{-1}}$ \\
% %                 \hfill $S=\rho g(\beta_s + \vartheta \beta_w)$ \\
% %           \begin{itemize}
% %             \item kompresibilita horniny $\beta_s = 0$ $\mathrm{[Pa^{-1}]}$
% %             \item kompresibilita vody $\beta_w = 4\cdot 10^{-10}$ $\mathrm{[Pa^{-1}]}$
% %             \item porozita $\vartheta = 0.007$ $\mathrm{[-]}$
% %           \end{itemize}
%     %
%         počáteční tlaková výška & $h_0$ & $300$ & $\unit{m}$ \\
%         Biotův-Willisův koeficient & $\alpha$ & $0.2$ & $\unit{-}$ \\
%         Youngův modul & $E$ & $6\cdot10^{10}$ & $\unit{Pa}$ \\
%         Poissonova číslo & $\nu$ & $0.2$ & $\unit{-}$ \\
%         \bottomrule
%     \end{tabular}
%     % TODO: zmenit popisek, diplnit citaci Rutquista (nyni je v textu), 
%     \caption{Parametry dopředného modelu.}
%     \label{tab:model_01_parameters}
% \end{table}

%Počáteční napětí v~horninovém masivu má rozpětí cca 20--35 MPa a poměr velikostí maximálního a minimálního horizontálního napětí se pohybuje okolo hodnoty 2.0, přičemž maximální horizontální napětí je obecně orientováno ve směru přibližně S-J až SV-JZ (\cite{Soucek2017Geotechnical}).

\begin{figure}[h]
    \centering
    \includegraphics[width=0.9\textwidth]{graphics/geometry_bc2.pdf}
    \caption{Schéma tunelu ve 2D řezu TSX experimentu.
    V~obrázku jsou znázorněny uvažované okrajové podmínky na stěně tunelu i na vnější hranici.
    Pro větší přehlednost nejsou celkové rozměry oblasti v~měřítku tunelu a pozic pozorovacích bodů.}
    \label{fig:model_01_geometry}
\end{figure}

Abychom zachytili změny v~EDZ, uvažujeme v poroelastickém modelu nelineární závislost mezi hydraulickou vodivostí a tenzorem napětí, viz rovnice \eqref{eqn:nonlinear_cond}. Konkrétně použijeme následující empirický vztah (\cite{Rutqvist2009Modeling}):

% \sy{[Standa: Nerozumím, jak je do poroelastického modelu zahrnut následující empirický vztah. Provádí se nejprve elastický výpočet, abychom z empirického vztahu získali vodivost potřebnou poroelastický model? Prosím o~podrobnější popis.]}
% - odkaz do kapitoly o~poroelasticite, tam upraveno od JS

\begin{equation} k = \frac{\rho g}{\mu} \left[\mathring{\permsc}_r + \Delta \mathring{\permsc}_m\exp(\mathring\beta\sigma_m)\right]\exp(\mathring\gamma\Delta\sigma_{VM}),
\label{eqn:nonlinear_cond_exp}
\end{equation}
kde
% TODO check translation
\begin{itemize}
    \item $\mathring \permsc_r=2\cdot10^{-21}\,\unit{m^2}$ je \emph{reziduální permeabilita při vysokém průměrném napětí}
    a fitovací konstanty: $\Delta \mathring{\permsc}_m=8\cdot10^{-17}\,\unit{m^2}$,  $\mathring\beta=4\cdot10^{-7}\,\unit{Pa^{-1}}$, $\mathring\gamma=3\cdot10^{-7}\,\unit{Pa^{-1}}$;
    \item $\sigma_m$ je průměr hlavních napětí, tj.
        \begin{equation}\label{eq:mean_stress} \sigma_m:=\frac13\tr\tn\sigma(\bf{u});\end{equation}
    \item $\sigma_{VM}$ je von Misesovo napětí:
        \begin{equation}\label{eq:von_mises} \sigma_{VM}:=\sqrt{\frac32\tn\sigma_d:\tn\sigma_d} 
            = \frac{1}{\sqrt{2}}\sqrt{\sqr{\sigma_1-\sigma_2}+\sqr{\sigma_2-\sigma_3}+\sqr{\sigma_1-\sigma_3}},\end{equation}
        kde $\tn\sigma_d:=\tn\sigma(\bf{u})-\frac13\big(\tr\tn\sigma(\bf{u})\big)\tn I$ je deviatorické napětí;
    \item $\Delta\sigma_{VM}$ je chápáno jako nezáporný rozdíl mezi vypočteným $\sigma_{VM}$ a
        \emph{kritickým deviatorickým napětím} $\sigma_{c}=55\,\unit{MPa}$, které odpovídá smykově indukované permeabilitě, viz \cite{Rutqvist2015Fractureda}. Platí tedy
         \begin{equation}\Delta\sigma_{VM} = \max\left(0, \sigma_{VM} - \sigma_{c}\right).\end{equation}
        Poslední člen ve vztahu \eqref{eqn:nonlinear_cond_exp} má význam jednoduchého modelu smykové dilatace. Přesnější popis tohoto efektu diskutujeme v~kapitole \ref{sec:nonlinear_fr_effects}.
 \end{itemize}


V~Tabulce \ref{tab:model_01_parameters} shrnujeme variabilní parametry dopředného modelu včetně uvažovaného apriorního rozdělení v~následné Bayesovské analýze. \uv{Střední hodnoty} těchto materiálových parametrů odpovídají doporučení od~\cite{Rutqvist2009Modeling}, které získali na základě kalibrace pro měřené pórové tlaky v~pozorovacím bodě V1.
\uv{Střední hodnotou} rozumíme hodnotu v uvedeném (srozumitelném) fyzikálním rozměru, která odpovídá střední hodnotě uvažovaného pravděpodobnostního rozdělení (např. logaritmus uvedené \uv{střední hodnoty} storativity je střední hodnota uvažovaného lognormálního rozdělení; pro beta rozdělení se rozumí střední hodnota ekvivalentního normálního rozdělení).
Zvolenou \uv{odchylku} pro apriorní lognormální rozdělení uvádíme v dekadickém logaritmu, v případě beta rozdělení je to odchylka ekvivalentního normálního rozdělení.

% Pro zachycení měřených porézních tlaků v~pozorovacím bodě V1 doporučili
% \cite{Rutqvist2009Modeling} hodnoty materiálových parametrů uvedených v~Tabulce \ref{tab:model_01_parameters}.    
% Hodnoty uvedené v~Tabulce \ref{tab:model_01_parameters} použijeme v~Bayesovské inverzní analýze pro a~apriorní odhad rozdělení parametrů.


\begin{table}[ht]
    \centering
    \begin{tabular}{p{3.0cm}crccc}
        % TODO: sloucit s Tab. 3
        % sloupecky "stredni hodnota", "odchylka", "rozdeleni" - popsat v textu transformaci na rozdeleni
        % logaritmus hodnoty je stredni hodnota lognormalniho rozdeleni
        % sigma v dekadickem logaritmu, samostatny sloupecek
        % biot: sigma bude odchylka ekvivalentniho normalniho rozdeleni
        \toprule
        \rot{parametr} & \rot{symbol} & \rot{\uv{střední hodnota}} & \rot{jednotka} & \rot{rozdělení} & \rot{\uv{odchylka}}\\
        \midrule
        Biotův-Willisův koeficient & $\alpha$ & $0.20$ & -- & beta & 0.025 \\
        storativita & $S$ & $7.29\cdot10^{-8}$ & $\unit{m^{-1}}$ & lognormální & 0.87\\
        Youngův modul & $E$ & $60$ & $\unit{GPa}$ & lognormální & 0.87 \\
        Poissonovo číslo & $\nu$ & $0.2$ & -- & pevné & -- \\[5pt]
        \multirow{4}{*}[-2pt]{hydraulická vodivost}
        %& $k$ & $5.89\cdot10^{-15}$ & $\unit{m\cdot s^{-1}}$ & lognormální & 3.0\\
        & $\mathring{\permsc}_r$ & $6.00\cdot10^{-22}$ & $\unit{m^2}$ & lognormální & 1.30 \\
        & $\mathring{\permsc}_m$ & $8.00\cdot10^{-17}$ & $\unit{m^2}$ & lognormální & 1.30 \\
        & $\mathring{\beta}$ & $4.00\cdot10^{-7}$ & $\unit{Pa^{-1}}$ & lognormální & 0.87 \\
        & $\mathring{\gamma}$ & $3.00\cdot10^{-7}$ & $\unit{Pa^{-1}}$ & lognormální & 0.87 \\
%                 \hfill $S=\rho g(\beta_s + \vartheta \beta_w)$ \\
%           \begin{itemize}
%             \item kompresibilita horniny $\beta_s = 0$ $\mathrm{[Pa^{-1}]}$
%             \item kompresibilita vody $\beta_w = 4\cdot 10^{-10}$ $\mathrm{[Pa^{-1}]}$
%             \item porozita $\vartheta = 0.007$ $\mathrm{[-]}$
%           \end{itemize}
    %
        % počáteční \\tlaková výška & $h_0$ & $300$ & $\unit{m}$ \\
        \bottomrule
    \end{tabular}
    % TODO: zmenit popisek, diplnit citaci Rutquista (nyni je v textu), 
    \caption{Vybrané parametry dopředného modelu pro Bayesovskou inverzi. Pro srovnání zpráva \cite{TZ303_2018}, Tab. 13, uvádí pro lokalitu Kraví hora Youngův modul v~rozmezí 36--50 GPa a Poissonovo číslo v~rozmezí 0.17--0.23.}
    \label{tab:model_01_parameters}
\end{table}

V modelu dále uvažujeme následující fixní hodnoty hustoty vody, gravitačního zrychlení a příčného řezu 2D oblasti:
\[\varrho_l=1000\,\unit{kg\cdot m^{-3}},\; g=9.81\,\unit{m\cdot s^{-2}} \mbox{ and } \delta=1\,\unit{m}.\]

%\subsection{Popis úlohy Bayesovské inverzní analýzy}

%\sy{[Standa: Trochu mě překvapuje lognormální rozdělení pro Youngův modul pružnosti nebo pro složku počátečního napětí. Tam bych zas až tak velké skoky jako u vodivosti nečekal.]}
% Jednak se jedná o prior a je vidět že v případě V1 dojde ke zůžení. Myslím také, že stejně jako pro vodivost je zde efekt závislosti tuhosti na měřítku. 
% Ovšem pokud bychom stanovili přesnější Youngův modul, tak se nám díky pozorované vazbě na storativitu a vodivost patrně upřesní i tyto související veličiny.
% Ohledně počátečního napětí, tam došlo pro V1 k ještě výraznějšímu zpřesnění prakticky na rozsah půl řádu.

% V Bayesovské inverzní analýze, popsané v~kapitole \ref{sec:Bayes},
% vezmeme jako dopředný problém $G$ výše popsaný poroelastický model počítaný ve Flow123d.
% Vybranou sadu parametrů $\vc u$, které budeme uvažovat jako nezávislé náhodné veličiny, jsme již uvedli v~Tabulce \ref{model_01_parameters},
% včetně volby apriorního rozdělení. Uvažujeme nelineární vztah pro hydraulickou 
% vodivost \eqref{eqn:nonlinear_cond_exp}, kde jako náhodné veličiny vystupují všechny 4 parametry.
% Všechna apriorní rozdělení uvažujeme lognormální, vyjma Biotova koeficientu, kde navrhujeme beta rozdělení.

% \begin{table}[!htb]
%     \centering
%     \begin{tabular}{lllc}
%     % parameter
%         \toprule
%         parametr &  apriorní rozdělení \\
%         \midrule
%         $\alpha$ & $\mathcal{B}(\alpha=50, \beta=200)$\\
%         $\log(S)$ & $\mathcal{N}(\mu=-16.4, \sigma=2.0)$\\
%         $\log(E)$ & $\mathcal{N}(\mu=24.8, \sigma=2.0)$\\
%         $\log(\sigma_x)$ & $\mathcal{N}(\mu=17.6, \sigma=2.0)$\\
%         $\log(\mathring{\permsc}_r)$ & $\mathcal{N}(\mu=-48.9, \sigma=3.0)$\\
%         $\log(\mathring{\permsc}_m)$ & $\mathcal{N}(\mu=-37.1, \sigma=3.0)$\\
%         $\log(\mathring{\beta})$ & $\mathcal{N}(\mu=-14.7, \sigma=2.0)$\\
%         $\log(\mathring{\gamma})$ & $\mathcal{N}(\mu=-15.0, \sigma=2.0)$\\
%         \bottomrule
%     \end{tabular}
%     \caption{Parametry a apriorní rozdělení.}
%     \label{tab:bayes_parameters}
% \end{table}

% Měřená data $\vc y$ pro všechny 4 pozorovací body jsme získali interpolací kubickým splinem z dat v~grafu na \figref{fig:TSX}. Časovou osu pozorování, tedy výstup funkce $G(\vc u)$, jsme diskretizovali s krokem 10 dní. Uvažujeme aditivní šum s~rozdělením v~podobě Gaussovského náhodného procesu s~korelační délkou 30 dní a~standardní odchylkou 100~m.


% Uvažujeme model s aditivním šumem $y=G\left(u\right)+\eta$. K úplné
% specifikaci Bayesovské úlohy potřebujeme dopředný model $G:\Real^{n}\rightarrow\Real^{m}$,
% apriorní rozdělení neznámých parametrů a rozdělení šumu.
% \begin{itemize}
% \item Dopředný model:
% \begin{itemize}
% \item počítán pomocí Flow123d (lineární poroelasticita, později s modifikacemi)
% \item vstupy (neznámé parametry): parametry vodivosti, storativita, Youngův
% modul, Biotův koeficient, vektor počátečního napětí, apod.
% \begin{itemize}
% \item V nejjednodušším případě homogenní materiál, ale lze rozdělit oblast
% na podoblasti s různými (konstantními) materiály, případně uvažovat
% materiál v podobě náhodného pole -> zvýšení počtu neznámých parametrů.
% \end{itemize}
% \item výstupy:
% \begin{itemize}
% \item The observed data were in-situ measurements of the pore pressure in
% two boreholes (named HGT1 and HGT2) during 365 days; specifically,
% four time series named HGT1-5, HGT1-4, HGT2-4, HGT2-3 taken from {[}Candler,
% Cournut, Dixon, 2002, Fig. 3.52{]}
% \item jedná se tedy o~4 časové řady
% \item provedena diskretizace v čase, tj. převedeno na vektor pozorovaných
% hodnot v následujících časech: \\
% 0.0 - 10.0 - 17.0 - 27.0 - 37.0 - 47.0 - 57.0 - 67.0 - 77.0 - 87.0
% - 97.0 - 100.0 - 120.0 - 140.0 - 160.0 - 180.0 - 200.0 - 220.0 - 240.0
% - 260.0 - 280.0 - 300.0 - 320.0 - 340.0 - 360.0 - 365.0
% \end{itemize}
% \end{itemize}
% \item Apriorní rozdělení:0
% \begin{itemize}
% \item parametry nezávislé
% \item dle typu buďto lognormální nebo beta rozdělení
% \begin{itemize}
% \item Beta použito pro Biotův koeficient, lognormální pro většinu ostatních
% parametrů
% \end{itemize}
% \end{itemize}
% \item Rozdělení šumu:
% \begin{itemize}
% \item ``Zde bylo zvoleno rozdělení v podobě Gaussovského náhodného procesu,
% což reflektuje očekávanou závislost mezi šumem blízkých údajů v časových
% řadách''
% \item Gaussovský proces, pro každou časovou řadu zvlášť, distance=30, corr\_length=100,
% autokorelační funkce:
% \begin{itemize}
% \item default - Ornstein-Uhlenbeck:\\
% exp(-distance/corr\_length)
% \item squared exponential:\\
% exp(-(distance{*}{*}2)/(2{*}corr\_length{*}{*}2))
% \end{itemize}
% \item (v našich úlohách měl tento typ šumu zanedbatelný vliv oproti nekorelovanému)
% \end{itemize}
% \end{itemize}

\subsection{Bayesovská inverze a numerické výsledky}

V Bayesovské inverzní analýze, popsané v~kapitole \ref{sec:Bayes},
vezmeme jako dopředný problém $G$ výše popsaný poroelastický model počítaný ve Flow123d.
Vybranou sadu parametrů $\vc u$, které budeme uvažovat jako nezávislé náhodné veličiny, jsme již uvedli v~Tabulce \ref{tab:model_01_parameters},
včetně volby apriorního rozdělení. Uvažujeme nelineární vztah pro hydraulickou 
vodivost \eqref{eqn:nonlinear_cond_exp}, kde jako náhodné veličiny vystupují všechny 4 parametry.
Všechna apriorní rozdělení uvažujeme lognormální, vyjma Biotova koeficientu, kde navrhujeme beta rozdělení.

Měřená data $\vc y$ pro všechny 4 pozorovací body jsme získali interpolací kubickým splinem z dat v~grafu na \figref{fig:TSX}. Časovou osu pozorování, tedy výstup funkce $G(\vc u)$, jsme diskretizovali s krokem 10 dní. Uvažujeme aditivní šum s~rozdělením v~podobě Gaussovského náhodného procesu s~korelační délkou 30 dní a~standardní odchylkou 100~m.

\figref{fig:bayes_observe_histograms} zachycuje histogram časového průběhu tlaku pro jednotlivé přijaté vzorky z Bayesovské inverze. Histogram je porovnán s měřeným průběhem tlaků (světle modrá) a vzorkem s nejlepším souhlasem s měřením (best fit - červená). Inverze je provedena zvlášť pro každé ze čtyř míst měření pórového tlaku. 

\begin{figure}[htb!]
    \centering
    \begin{subfigure}[t]{0.48\textwidth}
      \includegraphics[width=\textwidth]{\figpath bayes_hist_G1.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
      \includegraphics[width=\textwidth]{\figpath bayes_hist_G2.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
      \includegraphics[width=\textwidth]{\figpath bayes_hist_G3.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \begin{subfigure}[t]{0.48\textwidth}
      \includegraphics[width=\textwidth]{\figpath bayes_hist_G4.pdf}
      %\caption{}
      %\label{}
    \end{subfigure}
    \caption{Histogramy pozorovacích dat pro inverze s jednotlivými vrty.
        Na pozadí histogramů můžeme srovnat měření (světle modrá), nejlepší fit (červená)
        a~interkvartilové rozpětí (tečkovaně).
        % \sy{[Pozn. od EH: chybné jednotky tlaku na svislých osách grafů.]}
        }
    \label{fig:bayes_observe_histograms}
\end{figure}

U blízkého bodu v~nadloží V1
dostáváme dobrou shodu nejlepšího fitu s~měřenými daty, podobně jako v~\cite{Rutqvist2009Modeling}.
Zároveň vidíme, že interkvartilové rozpětí dobře pokrývá měřená data.
U vzdálenějšího bodu V2 vidíme relativně dobrý nejlepší fit, ale již pozorujeme, že měřená data
nejsou zcela uvnitř interkvartilového rozpětí. Histogram také nevykazuje pokles tlaku ke konci pozorovaného období. 

Podobně jako v~citovaném článku, i zde je problém zachytit průběhy tlaků v~horizontálním vrtu.
Nejlepší fit v~případě blízkého pozorovacího bodu H1 vykazuje dobrou shodu s daty, nicméně
v~histogramu takový trend nepozorujeme. Také měřená data nejsou ani v~mezích 5-95 percentilů.
Zřejmě zde také nedochází k~takovému překmitu do záporných tlaků, jako jsme viděli v~\figref{fig:TSX}.
Trend relativně rychlého návratu tlaku na původní hodnotu ve vzdálenějším bodě ukazuje na nějaký skrytý jev,
který dopředný model vůbec nepostihuje. Může jít např. o~neznámé puklinové propojení v~blízkosti měřicí sondy, kvůli kterému dojde k~nárůstu tlaku dříve, než by se očekávalo vzhledem k~malým hydraulickým vodivostem dále od tunelu.

% 836, 2749, 0, 3585
% 3871, 3027, 85530, 92428
% 1912, 1453, 100109, 103474

\begin{table}[!htb]
    \centering
    \begin{tabular}{lrrr}
        \toprule
        fáze &  \centering akceptované & zamítnuté & před-zamítnuté \\
        \midrule
        MH & 836 & 2749 & 0\\
        DAMH + (surrogate update) & 3871 & 3027 & 85530\\
        DAMH & 1912 & 1453 & 100109 \\
        \bottomrule
    \end{tabular}
    \caption{Fáze Bayesovské inverze a počty vzorků.}
    \label{tab:bayes_nsamples}
\end{table}

V Tabulce \ref{tab:bayes_nsamples} uvádíme počty vzorků počítaných v~jednotlivých fázích
Bayesovské inverze pro pozorovací bod V1. V~první fázi se vyvíjí zástupný model, ve druhé fázi se již používá pro
před-zamítnutí vzorků a dále se vyvíjí, v~poslední fázi se zástupný model pouze používá a zůstává neměnný.
Součet akceptovaných a zamítnutých vzorků představuje počet běhů poroelastického modelu ve Flow123d.
Pokud sečteme počty vzorků ve všech fázích dohromady, činí počet vyhodnocení zástupného modelu přes 96~\% všech vzorků. Celkový čas simulace při 14 paralelně počítaných řetězcích trval na clusteru Metacentra okolo 8~h.

\begin{figure}[htb!]
    \centering
    \includegraphics[width=\textwidth]{\figpath bayes_histograms2_G1G3_stor.pdf}
    \caption{Histogramy vybraných dvojic parametrů a srovnání blízkých pozorovacích vrtů V1 a H1.
    Na diagonále porovnání apriorního (oranžově) a posteriorního rozdělení (červeně).
    Mimo diagonálu histogramy ve 2D řezu parametrickým podprostorem, červeným bodem označen nejlepší fit.
    }
    \label{fig:bayes_histograms_stor}
\end{figure}

Na \figref{fig:bayes_histograms_stor} jsou zobrazeny histogramy ve 2D řezech parametrickým prostorem, tedy
vždy pro dvojici vybraných parametrů. Můžeme zde porovnat výsledky mezi blízkými pozorovacími body V1 a H1.
Vidíme zde zřejmé korelace mezi parametry v~nadloží tunelu -- nepřímou závislost Youngova modulu a storativity (korelační koef. 0.43)
a jedním z parametrů vodivosti \texttt{perm\_kr} (tj. $\mathring{k}_r$ v~\eqref{eqn:nonlinear_cond_exp}, korelační koef. 0.60), který charakterizuje neporušený masiv, dále přímou závislost mezi storativitou 
a tímto parametrem vodivosti. Naopak ve vrtu na boku tunelu žádné takové korelace nepozorujeme.
Posteriorní rozdělení parametru vodivosti v~případě bodu V1 je výrazně užší než zvolený prior,
máme tedy přesnější odhad (to platí i pro ostatní parametry).
V~případě H1 k~výraznému zpřesnění prioru nedochází v~žádném z vybraných parametrů.
Na \figref{fig:bayes_observe_histograms} jsme pozorovali rozumný nejlepší fit, tudíž model je schopen zachytit průběh tlaku v~H1.
Zřejmě jsme ale nastavili příliš vysokou úroveň šumu, a proto posterior vychází tak široký.
% pozorované hodnoty zřejmě nejsou dostatečně citlivé na změnu parametrů a dopředný model zřejmě neumí dobře zachytit
% měřená data, viz již histogram na \figref{fig:bayes_observe_histograms}.
% PE: upraveno
% \jb{JB: Není jasné co je myšleno  "v prvním a druhém případě" tam kde je dobrý best fit a špatný histogram je model dobrý, ale asi máme nastavený příliš velký šum.}
% \jb{ v celé kapitole nemluvit o pozorovacích vretech ale pozorovacích bodech, zavést vlstní značení, které si je člověk schopen zapamatovat např. V1, V2, H1, H2}

\begin{figure}[htb!]
    \centering
    \includegraphics[width=0.75\textwidth]{\figpath bayes_histograms2_G1G3_stress.pdf}
    \caption{Histogramy vybraných dvojic parametrů a srovnání blízkých pozorovacích vrtů V1 a H1.
    Na diagonále porovnání apriorního (oranžově) a posteriorního rozdělení (červeně).
    Mimo diagonálu histogramy ve 2D řezu parametrickým podprostorem, červeným bodem označen nejlepší fit.
    }
    \label{fig:bayes_histograms_stress}
\end{figure}

Na \figref{fig:bayes_histograms_stress} můžeme porovnat další dva parametry.
Počáteční napětí ve směru $x$ je v~případě pozorovacího bodu V1 velmi blízko
hodnotě určené v~\cite{Rutqvist2009Modeling}. V~případě bočního pozorovacího bodu H1 posterior
ukazuje vyšší pravděpodobnost nižšího počátečního napětí.
Parametr vodivosti \texttt{perm\_beta} (tj. $\mathring\beta$ v~\eqref{eqn:nonlinear_cond_exp}) je v~případě bodu V1 spíše vyšší,
naopak v~případě bodu H1 je v grafu hustota pravděpodobnosti posunuta k nižším hodnotám.
To naznačuje další možnou heterogenitu, která není dopředným modelem zachycena.
% PE: upraveno
% \jb{Není jasné co je parametr $\beta$, použít raději i označení z obrázku "perm kr". Není jasné co znamená posunul se na druhou stranu. }

Rozptyl návrhové hustoty $q(x|k)$ pro kroky náhodné procházky bylo nutné volit dostatečně malý,
% PE: upresneno
%\jb{ velké co? Odpovídal by velká rozptyl, ale dále se tvrdí že je odchylka volena standardně??}
jinak docházelo k~tomu, že se algoritmus v~určitém bodě zastavil a všechny další vzorky zamítal.
Pro jednotlivé fáze jsme nakonec nastavili standardní odchylku rozdělení návrhové hustoty 0.1, 0.2 a 0.3.

Otevřenou otázkou zůstává \uv{rozumné} nastavení šumu, kterým kompenzujeme jednak chyby měření a jednak nedostatky dopředného modelu. Prozatím byla pro všechny pozorovací vrty volena stejná úroveň šumu, a to konstantní v~průběhu času. Zde by bylo potřeba provést více experimentů pro porovnání a~šum různě upravit pro jednotlivá pozorování a případně i změnit v~časové ose.


Dosavadní výsledky Bayesovské inverzní analýzy jsou v dobré shodě s~článkem \cite{Rutqvist2009Modeling} v případě fitování pozorovacích bodů v nadloží. 
V~případě pozorovacích bodů v bočním vrtu je nalezený "best fit" pro měření H1 výrazně lepší, avšak obecně hodnoty parametrů z posteriorního rozdělení nevykazují lepší shodu s měřením. Je zřejmé, že dopředný model je při vhodné volbě parametrů ve velmi dobré shodě s měřením a šum měření je tak zbytečně nadhodnocen. Pro zvládnutí nižšího šumu však bude třeba dále vylepšit algoritmus Bayesovské inverze o~adaptivní volbu návrhové hustoty, na kterou je algoritmus poměrně citlivý. Pro zachycení všech čtyř měření jednou množinou parametrů bude nutné uvažovat jejich heterogenitu. Hlavním přínosem Bayesovské inverze  je získání pravděpodobnostního rozdělení parametrů (především vodivosti a porozity), které nám umožňují dále propagovat stochastický přístup a zahrnout vliv EDZ v transportním modelu.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Závěrečné shrnutí}
\label{sec:zaver}

V předložené metodice jsme navrhli indikátor bezpečnosti vhodný zejména pro relativní porovnání bezpečnosti jednotlivých pozic ÚOS. Popsali jsme také výpočet indikátoru pomocí vhodného transportního modelu. Pro lepší interpretaci by byl vhodnější inikátor bezpečnosti založený na výsledné dávce, avšak jeho numerická realizace by vyžadovala vhodné aproximace dalších transportních procesů. 

Demonstrovali jsme realizaci stochastického výpočtu s náhodnými malými puklinami a aproximaci hustoty pravděpodobnosti indikátoru bezpečnosti. Výpočet byl proveden celkem pro 12 variant: dvě varianty uspořádání deterministických puklin, tři polohy zdroje kontaminace a varianta s a bez zvýšené vodivosti na EDZ. Z prezentovaných výsledků vyplývá velká citlivost indikátoru bezpečnosti  na vodivost v EDZ. Citlivost na polohu zdroje kontaminace byla pozorována pouze pro případ bez EDZ, obecněji tedy pro dostatečně nízkou hydraulickou vodivost na EDZ.

Na základě porovnání uvedených variant jsme vytvořili pracovní 
hypotézu pro další výzkum, že s rostoucí hydraulickou vodivostí 
EDZ roste počet puklin zapojených do transportu kontaminace
směrem k malým puklinám s malou vodivostí.
Tím roste i pravděpodobnost vytvoření preferenční cesty.



Použité stochastické metody je možno snadno modifikovat také pro získání korelací mezi hodnotami indikátoru resp. koncentrace v různých časech nebo korelací mezi indikátory pro různé úložné pozice. Bylo by tak možné získat kompletní stochastický popis výstupní koncentrace, která pak může být vstupem do odděleného modelu pro transport v geologické bariéře a navazujícího modelu šíření v~biosféře. 

Pro stochastický popis klíčových parametrů transportního modelu: hydraulické vodivosti a porozity na EDZ
jsme použili metodu Bayesovské inverze. Tuto stochastickou inverzní metodu jsme konkrétně aplikovali na kontinuální měření pórových tlaků v průběhu a po dokončení ražby. Za tímto účelem jsme vyvinuli hydro-mechanický model popisující změny pórového tlaku, hydraulické vodivosti a porozity 
v důsledku změn napjatosti horniny způsobenou ražbou. Demonstrovali jsme zvládnutí Bayesovské inverze 
i v kombinaci s výpočetně relativně náročným dopředným modelem. Toto bylo prakticky umožněno zejména díky implementaci Metropolis-Hastings metody se zpožděným přijetím v kombinaci se zástupnými modely. 

Použití metodiky je demonstrováno pouze pro variantu horizontálního ukládání s mechanickou ražbou úložných horizontálních vrtů. Většinu modelů lze poměrně snadno upravit pro charkterizaci EDZ v okolí přístupových chodeb s úložnými vertikálními vrty. Velmi nevýhodný profil chodeb by mohl vyžadovat použití plastických modelů (viz kapitola 
\ref{sec:plasticita} dodatku). Numerickými modely je však těžko postihnutelný vliv trhacích prací při použití klasické ražby. 
Z hlediska EDZ nelze systém vertikálního ukládání doporučit, avšak pro kvantifikaci vlivu systému ukládání na indikátor bezpečnosti by bylo třeba provést podrobnější studii, která je s využitím metodiky a zejména vyvinutého software jistě realizovatelná. 

Navržená metodika je prakticky realizovatelná pomocí komplexního výpočetního software, který je spolu s metodikou také výsledkem výzkumného projektu. V kontextu modelování je relativně unikátní sestavení parametrického geometrického modelu. Při změně parametrů 
geometrie není nutné ručně zpracovávat nový výpočetní model. Tento přistup tak umožňuje snadné porovnávání různých variant.

Software využívá řadu různých výpočetních nástrojů a knihoven a umožňuje využití paralelizace zejména pro stochastické výpočty. Zvládnutí komplexního výpočtu tak otevírá možnost sestavení celkového bezpečnostního modelu s~využitím otevřeného softwaru. Absence omezení vyplývajících z komerčních licencí a uzavřených datových formátů umožňuje poměrně snadnou rozšiřitelnost a propojitelnost různých výpočetních komponent.


\subsection{Směry dalšího vývoje}
Předkládanou metodiku je třeba chápat jako koncept, který je pro praktické nasazení nutné dále rozpracovat a validovat. 
V následujících odstavcích stručně shrneme oblasti, kterým je dle našeho mínění třeba věnovat prioritní pozornost. 

Je třeba dokončit aplikaci homogenizačních algoritmů pro složité geometrie a demonstrovat zrychlení stochastického výpočtu pomocí víceúrovňové metody Monte Carlo. Pro zlepšení robustnosti bude vhodné automatizovat volbu různých numerických parametrů, např. volbu časového kroku a volbu ukončovacích tolerancí vnitřních lineárních řešičů.  

Demonstrace metodiky byla realizována s využitím dílčích dat z různých lokalit, neboť pro žádnou dostupnou lokalitu nebyla k dispozici všechna potřebná měření. V příloze \ref{sec:Bukov} je proto představen návrh experimentální ražby s řadou doprovodných měření. Na základě takového experimentu bude možné validovat metodiku na reálných datech. Ačkoliv metoda kontinuálních měření pórových tlaků umožňuje poměrně přesně stanovit řadu parametrů horniny, její aplikace je příliš náročná pro plošné použití v rámci budování úložiště. Je proto potřeba pomocí porovnání s neinvazivními geofyzikálními metodami zjistit, zda lze variabilitu parametrů v rámci důlního díla charakterizovat pomocí geofyzikálních metod.

Bude vhodné doplnit použité dopředné modely. Transportní model již obsahuje implementaci sorpcí a radioaktivního rozpadu, je však třeba mít k dispozici parametry těchto modelů a stanovení množiny studovaných radionuklidů. Upřesnění je potřebné také pro zdrojový člen. Pro hydro-mechanický model bude vhodné zahrnutí plasticity pro popis nevratných mechanických jevů na mikro puklinách. V rámci projektu byl rozpracován samostatný plastický model, viz část \ref{sec:plasticita} přílohy. Mikroskopický model 
mechanických procesů na puklinách a~zvládnutí rozsáhlých úloh s kontakty a třením by měly být možné s využitím postupů popsaných v části \ref{sec:advaced_mech} přílohy. Pomocí víceškálového přístupu
by bylo možné využít přesnější konstitutivní vztah pro závislost tenzoru hydraulické vodivosti na tenzoru deformací, viz sekci \ref{sec:hm_homogenizace}. 

Plný potenciál Bayesovské inverze nebyl doposud využit zejména kvůli relativní citlivosti na 
volbu rozptylu návrhové hustoty. Ukazuje se, že pro plné využití relativně přesných měření je vhodné měnit adaptivně návrhovou hustotu v průběhu algoritmu. Dále bude třeba pomocí vhodné parametrizace zvládnout inverzi pro prostorově heterogenní parametry.  Dosavadní experimenty v tomto směru pouze ukázaly, že pro získání dobré shody ve všech bodech měření pórového tlaku nestačí uvažovat pouze heterogenitu hydraulické vodivosti, ale je nutné zahrnout také heterogenitu dalších proměnných.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\appendix

\section{Možnosti rozšíření metodiky}
\label{sec:rozsireni}

V této kapitole diskutujeme další možnosti rozšíření metodiky, které zahrnují: a) numerickou predikci EDZ, b) užití pokročilých numerických metod řešení poroelastických modelů s puklinami i bez nich a c) návrh homogenizačních technik pro určení anizotropie tenzoru hydraulické vodivosti na základě informace o~puklinové síti. 

\subsection{Predikce EDZ pomocí pružných a pružně-plastických modelů} 
\label{sec:plasticita}
Matematické modelování EDZ může být založeno na mechanice kontinua, lomové mechanice nebo na víceškálových přístupech. My se zaměříme na kontinuální modely, konkrétně na pružné a pružně-plastické modely. Od těchto modelů obvykle očekáváme robustnost vůči časové a prostorové diskretizaci, což by mělo být jejich přirozenou výhodou oproti plně diskrétním modelům. Pružné modely jsou obvykle kombinovány s kritériem porušení, díky kterému detekujeme zóny s vysokou koncentrací napětí. Tento přístup je nejjednodušší a je vhodný pro rozsáhlé 3D geometrie. V~případě elasto-plasticity je kritérium poškození přímočaře zahrnuto do modelu a přípustná pole napětí musí toto kritérium splňovat. Pružně-plastické modely mohou být obohaceny o~vnitřní proměnné představující změkčení, zpevnění nebo další vnitřní pomocné proměnné. V~článku \cite{perras_predicting_2016} a jeho citacích byly klasifikovány různé typy zón poškození určených na základě pružně-plastických modelů, tzv. DISL přístup. 

\subsubsection{Numerické schéma pružných a pružně-perfektně plastických modelů}
\label{sec_plasticity}

\begin{figure}[h]
	\begin{center}
	\begin{picture}(120,120)
	%ctverec
	\drawline(16,16)(116,16)(116,116)(16,116)(16,16)
	
	%srafovani
	{\linethickness{0.1mm}
		\drawline(16,106)(26,116)
		\drawline(16,96)(36,116)
		\drawline(16,86)(46,116)
		\drawline(16,76)(56,116)
		\drawline(16,66)(66,116)
		%\drawline(16,56)(76,116)%nahrada
		\drawline(67,107)(76,116)
		\drawline(16,56)(55,95)
		\drawline(16,46)(86,116)
		%\drawline(16,36)(96,116)%nahrada
		\drawline(16,36)(45,65)
		\drawline(60,80)(96,116)
		%\drawline(16,26)(106,116)%nahrada
		\drawline(16,26)(45,55)
		\drawline(70,80)(106,116)
		%\drawline(16,16)(116,116)%nahrada
		\drawline(16,16)(50,50)
		\drawline(80,80)(116,116)
		%\drawline(26,16)(116,106)%nahrada
		\drawline(26,16)(60,50)
		\drawline(85,75)(116,106)
		%\drawline(36,16)(116,96)%nahrada
		\drawline(36,16)(70,50)
		\drawline(95,75)(116,96)
		%\drawline(46,16)(116,86)%nahrada
		\drawline(46,16)(88,58)
		\drawline(105,75)(116,86)
		%\drawline(56,16)(116,76)%nahrada
		\drawline(56,16)(64,24)
		\drawline(75,35)(98,58)
		\drawline(110,70)(116,76)
		\drawline(66,16)(116,66)
		\drawline(76,16)(116,56)
		\drawline(86,16)(116,46)
		\drawline(96,16)(116,36)
		\drawline(106,16)(116,26)}
	
	%trojuhelnicky a kolecka dole
	\multiput(0,0)(29,0){4}{\drawline(18,4)(28,4)(23,12)(18,4)}
	\multiput(0,0)(29,0){4}{\put(20,0){\circle{4}}}
	\multiput(0,0)(29,0){4}{\put(26,0){\circle{4}}}
	
	%trojuhelnicky a kolecka nahore
	\multiput(0,0)(29,0){4}{\drawline(18,128)(28,128)(23,120)(18,128)}
	\multiput(0,0)(29,0){4}{\put(20,132){\circle{4}}}
	\multiput(0,0)(29,0){4}{\put(26,132){\circle{4}}}
	
	%trojuhelnicky a kolecka vlevo
	\multiput(0,0)(0,29){4}{\drawline(4,18)(12,23)(4,28)(4,18)}
	\multiput(0,0)(0,29){4}{\put(0,20){\circle{4}}}
	\multiput(0,0)(0,29){4}{\put(0,26){\circle{4}}}
	
	%trojuhelnicky a kolecka vpravo
	\multiput(0,0)(0,29){4}{\drawline(128,18)(120,23)(128,28)(128,18)}
	\multiput(0,0)(0,29){4}{\put(132,20){\circle{4}}}
	\multiput(0,0)(0,29){4}{\put(132,26){\circle{4}}}
	
	%sipka,oval
	\put(66,66){\oval(32,22)}
	\put(99,65){\vector(-1,0){15}}
	
	%popis
	\put(69,30){\makebox(0,0){$\Omega$}}
	\put(60,99){\makebox(0,0){$\mbf\sigma_0$}}
	%	\put(65,65){\makebox(0,0){$\Omega_r$}}
	\put(105,65){\makebox(0,0){$\Gamma$}}
	
\end{picture}
\caption{2D geometrie horninového masivu v okolí tunelu.}
\label{fig_scheme}
\end{center}
\end{figure}

Uvažujme zjednodušenou 2D geometrii, která představuje schéma horninového masivu v~okolí tunelu, viz \figref{fig_scheme}. Čtvercová oblast bude značena symbolem $\hat\Omega$, zatímco její vyšrafovaná podoblast  $\Omega$. Na vnější hranici $\partial\hat\Omega$ (daleko od tunelu) předepisujeme nulové normálové posuny a na vnitřní hranici $\Gamma$ uvažujeme nulové normálové napětí. Tyto okrajové podmínky zapisujeme ve tvaru $\mbf u\cdot\mbf n=0 $ na $\partial \hat\Omega$ a $\mbf\sigma\mbf n=\mbf0$ na $\Gamma$, kde $\mbf u$, $\mbf\sigma$ a $\mbf n$ označují respektive pole posunů, pole napětí a vnější jednotkový normálový vektor na $\Omega$. Dále předepisujeme počáteční pole napětí $\mbf\sigma_0$ definované v~$\Omega$. Pro zjednodušení simulujeme ražbu tunelu pomocí historie zatížení zapsané ve tvaru $t\mbf\sigma_0/t_{max}$, kde parametr $t\in[0,t_{max}]$ reprezentuje zkoumaný časový úsek. Dalšími složkami pružných a~pružně-plastických modelů je tenzor malých deformací
$$\mbf\varepsilon:=\mbf\varepsilon(\mbf u)= \frac{1}{2}(\nabla \mbf u + (\nabla\mbf u)^T)$$
a elastický tenzor čtvrtého řádu $\tn C$ splňující
\begin{eqnarray*}
	\tn{C} \mbf \varepsilon &=& \frac{E}{1+\nu}\left\{\frac{\nu}{1-2\nu} (\text{tr}\ \mbf \varepsilon)\tn I + \mbf \varepsilon \right\}, \label{C2}\\
	\tn{C}^{-1} \mbf \sigma &=& - \frac{\nu}{E} (\text{tr}\ \mbf \sigma)\tn I + \frac{1+\nu}{E} \mbf \sigma,
	\label{C_inv2}
\end{eqnarray*}
kde $E>0$, $\nu\in(0,1/2)$ značí respektive Youngův modul a Poissonovo číslo, $\tn I$ je jednotkový tenzor druhého řádu a $\text{tr}\ \mbf \varepsilon=\mbf \varepsilon:\tn I=\varepsilon_{11}+\varepsilon_{22}+\varepsilon_{33}$ značí stopu tenzoru $\mbf\varepsilon$.

V případě lineární pružnosti máme následující konstitutivní (Hookův) zákon mezi tenzory napětí a malých deformací:
$$\mbf \sigma=\tn{C}\mbf\varepsilon + \mbf\sigma_0\quad\mbox{or}\quad \mbf \sigma=\tn{C}[\mbf\varepsilon + \mbf\varepsilon_0],\;\;\mbf\varepsilon_0
=\tn C^{-1}\mbf\sigma_0.$$

{\sloppy
Nyní uvedeme pružně-perfektně plastický konstitutivní model, který je oproti pružnému modelu časově závislý. Nechť $\mbf \varepsilon^e$ a $\mbf \varepsilon^p$ značí elastickou a plastickou část tenzoru malých deformací a $\lambda$ je takzvaný plastický multiplikátor. Dále potřebujeme definovat funkci poškození $f:=f(\mbf\sigma)$ a plastický potenciál $g:=g(\mbf\sigma)$ a předpokládáme, že obě funkce jsou konvexní. Potom konstitutivní problém formulujeme následovně: pro každé $t\in(0,t_{max})$, hledáme $\mbf \sigma:=\mbf \sigma(t)$, $\mbf \varepsilon:=\mbf \varepsilon(t)$, $\mbf \varepsilon^e:=\mbf \varepsilon^e(t)$, $\mbf \varepsilon^p:=\mbf \varepsilon^p(t)$, $\lambda:=\lambda(t)$ takové, že
}
\begin{itemize}
	\item $\mbf\varepsilon = \mbf \varepsilon^e + \mbf \varepsilon^p$, $\mbf\sigma = \tn{C}(\mbf\varepsilon^e + t\mbf\varepsilon_0/t_{max})$,
	\item $\dot{\mbf\varepsilon}^p \in \dot{\lambda}\partial g(\mbf\sigma)$, $\mbf\varepsilon^p(0)=0$,
	\item $\dot{\lambda} \geq 0$, $\dot{\lambda} f(\mbf\sigma) = 0$, $f(\mbf\sigma) \leq  0$.
\end{itemize}
Zde tečka nad symboly značí časovou derivaci a $\partial g(\mbf\sigma)$ znamená subdiferenciál funkce $g$ v~$\mbf\sigma$. Je důležité zmínit, že subdiferenciální definice plastického tečení není příliš obvyklá v~inženýrské praxi, ve které se běžně používá derivace funkce $g$: 
$$\dot{\mbf\varepsilon}^p = \dot{\lambda}\frac{\partial g(\mbf\sigma)}{\partial\mbf\sigma},$$
navzdory skutečnosti, že potenciál $g$ není všude diferencovatelný. Jedním z cílů našeho výzkumu bylo navázat na články \cite{sysala_subdifferential-based_2016, sysala_subdifferential-based_2017} a ukázat, že explicitní znalost tvaru subdiferenciálu funkce $g$ usnadňuje analýzu a numerické řešení pružně-plastického konstitutivního problému.

Pružně-plastický konstitutivní problém je většinou diskretizován v~čase pomocí implicitní Eulerovy metody. Uvažujeme-li následující dělení časového intervalu 
$0=t_0<t_1<\ldots<t_{k-1}<t_k<\ldots<t_{max}$, potom má diskretizovaný konstitutivní problém následující tvar: k~daným hodnotám $\mbf\varepsilon_0$, $\mbf\varepsilon_k$ a $\mbf \varepsilon_{k-1}^p$ hledáme $\mbf\sigma_k$ a $\mbf\varepsilon_{k}^p$ takové, že
$$\mbf \sigma_k = T(\mbf\varepsilon_k- \mbf \varepsilon_{k-1}^p+t_k\mbf \varepsilon_0/t_{max}),\quad \mbf \varepsilon_{k}^p = \mbf\varepsilon_k  + t_k\mbf \varepsilon_0/t_{max} - \tn{C}^{-1} \mbf \sigma_k.$$
Zde $T$ je neznámá tenzorová funkce, kterou potřebujeme zkonstruovat. Obecně je tato funkce daná pouze implicitním předpisem, a proto její konstrukce musí být provedena numericky, pomocí algoritmu zahrnujícího elastickou predikci a plastickou korekci, viz \cite{neto_computational_2008, clausen_exact_2008, sysala_subdifferential-based_2016, sysala_subdifferential-based_2017}. V~rámci elastické predikce ověřujeme, zda-li je tak zvané testovací napětí $\mbf \sigma_k^{tr} = \tn{C}(\mbf\varepsilon_k- \mbf \varepsilon_{k-1}^p+t_k\mbf \varepsilon_0/t_{max})$ plasticky přípustné, tj. jestli splňuje podmínku $f(\mbf \sigma_k^{tr})\leq0$ reprezentující plastické kritérium. Pokud je tato podmínka splněna, potom platí $\mbf \sigma_k=\mbf \sigma_k^{tr}$ a $\mbf \varepsilon_{k}^p =\mbf \varepsilon_{k-1}^p$. V~opačném případě aplikujeme plastickou korekci, při které musí být splněna podmínka $f(\mbf \sigma_k)=0$. To znamená, že potřebujeme navrátit predikované napětí na plochu plasticitu. Plastické korekci se proto často říká návratový algoritmus (anglicky return-mapping algorithm). Z matematického úhlu pohledu plastickou korekci interpretujeme jako zobecněnou projekci na konvexní množinu, díky čemuž odvozujeme důležité vlastnosti funkce $T$. Například platí, že tato implicitní funkce je semihladká a její hodnoty jsou jednoznačně definované, viz \cite{sysala_subdifferential-based_2016, sysala_subdifferential-based_2017}.

Užitím funkce $T$ můžeme celkový pružně-plastický problém formulovat vůči neznámému poli posunů v~následujícím tvaru:
$$\mbox{hledáme }\mbf u_k \in V:\;\;  \int_{\Omega} T(\varepsilon(\mbf u_k)- \mbf \varepsilon_{k-1}^p+t_k\mbf \varepsilon_0/t_{max}): \mbf\varepsilon(\mbf v)\, dx
= 0 \quad \forall\mbf v \in V,$$
kde
$$V=\{\mbf v\in H^1(\Omega;\Real^d)\ |\;\; \mbf v \cdot \mbf n =0 \ \text{on}\ \partial \hat\Omega \}$$
je prostor přípustných polí posunů. Po prostorové diskretizaci přicházíme k~řešení systému nelineárních rovnic. Tento systém je pak řešen nehladkou variantou Newtonovy metody. Abychom tuto metodu mohli použít, potřebujeme navíc zkonstruovat zobecněnou derivaci funkce $T$. Její konstrukci představíme v~další sekci pro uvažované plastické modely.


\subsubsection{Mohr-Coulombův a Hoek-Brownův konstitutivní model}
\label{sec_HB}

Mohr-Coulombovy a Hoek-Brownovy plastické modely jsou běžně používané v~inženýrské geotechnice. Funkce $f$ a $g$ jsou pro tyto modely definovány pomocí hlavních napětí, které jsou určeny na základě následujícího spektrálního rozkladu tenzoru napětí:
\begin{equation*}
	\mbf\sigma=\sum_{i=1}^3 \sigma_i\mathbf e_i\otimes\mathbf e_i,\quad \sigma_1\geq \sigma_2\geq \sigma_3.
	\label{spectrum}
\end{equation*}
Zde $\sigma_i\in\Real$  a $\mathbf e_i\in\Real^3$, $i=1,2,3$, značí vlastní čísla (hlavní napětí) a vlastní vektory. Předpokládáme, že hlavní napětí jsou uspořádány tak, že platí $\sigma_1\geq\sigma_2\geq \sigma_3$. V~dalším budeme pracovat s běžnou mechanickou znaménkovou konvencí, která předpokládá kladné hodnoty pro tahy, ačkoliv v~geotechnice se používá opačná znaménková konvence. 

Mohr-Coulombův model je definován pomocí funkcí
\begin{eqnarray*}
	f(\mbf\sigma)&=&(1+\sin\phi)\sigma_1-(1-\sin\phi)\sigma_3-2c\cos\phi,\\
	g(\mbf\sigma)&=&(1+\sin\psi)\sigma_1-(1-\sin\psi)\sigma_3,
	\label{potential_function}
\end{eqnarray*}
kde $c>0$, $\phi\in(0,\pi/2)$ a $\psi\in(0,\pi/2)$ jsou dané materiálové parametry značící kohezi, úhel tření a úhel dilatace. Pro většinu materiálů platí $\psi\leq\phi$. 

Hoek-Brownův model je definován pomocí funkcí
\begin{eqnarray*}
	f(\mbf\sigma)&=&\sigma_1-\sigma_{ci}\Big(s-m_b\frac{\sigma_1}{\sigma_{ci}}\Big)^a-\sigma_3, \\
	g(\mbf\sigma)&=&\sigma_1-\sigma_{ci}\Big(s_g-m_g\frac{\sigma_1}{\sigma_{ci}}\Big)^{a_g}-\sigma_3,
	%\label{potential_function}
\end{eqnarray*}
kde $\sigma_{ci}, s, s_g, m_b, m_g>0$ a $a,a_g\in(0,1)$ jsou dané materiálové parametry. \sy{Vidíme, že funkce $f$ je definována tzv. vrcholovými parametry $s, m_b, a$, zatímco funkce $g$ je určena pomocí tzv. reziduálních parametrů $s_g, m_g, a_g$. Pro křehké horniny jsou obvykle hodnoty reziduálních parametrů mnohem nižší než hodnoty vrcholových parametrů. Aby byl model dobře definovaný, odvodili jsme předpoklad $s_g/m_g\geq s/m$, který v~běžné geotechnické literatuře není uveden. V~limitním případě $a=a_g=1$ přecházíme od Hoek-Brownova modelu k~Mohr-Coulombovu modelu. Dále je důležité zmínit, že k definici materiálových parametrů pro horninový masiv se používají empirické vztahy, které obsahují běžné materiálové parametry pro intaktní horniny a dva přídavné parametry (indexy): geologický pevnostní index $GSI$ reprezentující strukturu okolního horninového masivu a index disturbance $D$ charakterizující způsob ražby tunelu či vrtu. Více detailů k Hoek-Brownovu modelu} lze najít v~\cite{clausen_exact_2008, hoek_hoekbrown_2019, hoek_hoek-brown_2002}. Stručně řečeno,   

Připomeňme, že přípustná pole napětí splňují podmínku $f(\mbf\sigma)\leq 0$. Množina všech přípustných napětí pro Mohr-Coulombův nebo Hoek-Brownův model je konvexní a zarovnaná podél hydrostatické osy (na které platí $\sigma_1=\sigma_2=\sigma_3$). Konkrétně, množina Mohr-Coulombových přípustných napětí má tvar šestibokého jehlanu (pyramidy) v~prostoru hlavních napětí, přičemž vrchol pyramidy je dán hydrostatickým napětím $\sigma_t=c/\tan\phi$. Pro Hoek-Brownův model, je šestiboká pyramida zakřivená a má vrchol v~bodě  $\sigma_t=s\sigma_{ci}/m_b$, viz \cite{clausen_exact_2008}. Dále si můžeme všimnout, že funkce $g$ lze pro oba modely zapsat následujícím jednotným způsobem:
$$g(\mbf\sigma)=\hat g_1(\sigma_1)-\hat g_3(\sigma_3),$$
kde $\hat g_1$ a $\hat g_3$ jsou rostoucí, konvexní a dvakrát diferencovatelné funkce na svém oboru hodnot.
Rozšíříme-li výsledky prezentované  v~\cite{sysala_subdifferential-based_2017}, je možné ukázat, že funkce $g$ je pro oba modely konvexní a její subdiferenciál má následující tvar:
\begin{eqnarray*}
	\partial g(\mbf\sigma)&=&\Big\{ \mbf\nu=\sum_{i=1}^3\nu_i\mathbf e_i\otimes\mathbf e_i\ \Big|\;\; (\mathbf e_1, \mathbf e_2, \mathbf e_3)\in V(\mbf\sigma); \; \nonumber\\
	&&\quad \hat g_1'(\sigma_1)\geq\nu_1\geq\nu_2\geq\nu_3\geq -\hat g'_3(\sigma_3);\;\;\sum_{i=1}^3\nu_i=\hat g_1'(\sigma_1)-\hat g'_3(\sigma_3);\nonumber\\
	&&\quad (\nu_1-\hat g_1'(\sigma_1))(\sigma_1-\sigma_2)=0;\; (\nu_3+\hat g'_3(\sigma_3))(\sigma_2-\sigma_3)=0 \Big\},
	\label{sub_g}
\end{eqnarray*}
kde $\hat g_1'$, $\hat g_3'$ značí derivace funkcí $\hat g_1$, $\hat g_3$ a

	
$$V(\mbf\sigma)=\{(\mathbf e_1,\mathbf e_2,\mathbf e_3)\in [\Real^3]^3\ |\; \mathbf e_i\cdot\mathbf e_j=\delta_{ij};\; \mbf\sigma\mathbf e_i=\sigma_i\mathbf e_i,\; i,j =1,2,3;\; \sigma_1\geq \sigma_2\geq \sigma_3\}$$
je prostor vlastních vektorů. Pokud platí $\sigma_1>\sigma_2>\sigma_3$, potom $\nu_1=\hat g_1'(\sigma_1)$, $\nu_2=0$ a $\nu_3=-\hat g'_3(\sigma_3)$, z čehož vyplývá, že funkce $g$ je diferencovatelná v~$\mbf\sigma$. V~opačném případě ale derivace funkce $g$ v~bodě $\mbf\sigma$ neexistuje a je zapotřebí pracovat pouze s jejím subdiferenciálem, přičemž hodnoty $\nu_1,\nu_2,\nu_3$ nejsou jednoznačně definovány.  

Připomeňme dále, že v~případě plastické korekce platí, že hledaný tenzor napětí $\mbf\sigma:=\mbf\sigma_k$ splňuje $f(\mbf\sigma)= 0$. V~takovém případě $\mbf\sigma$ leží na povrchu Mohr-Coulombovy nebo Hoek-Brownovy pyramidy. Vzhledem k~předepsanému uspořádání $\sigma_1\geq \sigma_2\geq \sigma_3$ můžeme povrch pyramidy rozdělit na čtyři části: hladkou část ($\sigma_1>\sigma_2>\sigma_3$), levou (zakřivenou) hranu ($\sigma_1=\sigma_2>\sigma_3$), pravou (zakřivenou) hranu ($\sigma_1>\sigma_2=\sigma_3$) a na vrchol ($\sigma_1=\sigma_2=\sigma_3=\sigma_t$). Tato terminologie byla zavedena v~\cite{neto_computational_2008}. Pro každý z těchto čtyř případů můžeme upřesnit tvar množiny $\partial g(\mbf\sigma)$ a následně i definici plastické korekce. Například, pokud nastává návrat na levou hranu, potom platí $\hat g_1'(\sigma_1)\geq\nu_1\geq\nu_2\geq0$, $\nu_1+\nu_2=\hat g_1'(\sigma_1)$,  $\nu_3=-\hat g'_3(\sigma_3)$. Tyto podmínky nejsou obvyklé v~inženýrské praxi, nicméně mohou zjednodušit konstrukci plastické korekce nebo odvodit, na jaké části povrchu pyramidy bude ležet hledané napětí.

V případě pružně-perfektně plastického Mohr-Coulombova modelu můžeme odvodit rozhodovací kritéria pro každý typ návratu na povrch pyramidy a následně najít explicitní tvar konstitutivní operátoru $T$, viz \cite{neto_computational_2008,sysala_subdifferential-based_2017}. Avšak pro Hoek-Brownův model je tenzorová funkce $T$ daná pouze implicitně. V~článku \cite{clausen_exact_2008} byl navržen následující algoritmus pro plastickou korekci. Nejprve testujeme, zda-li nenastává návrat do vrcholu Hoek-Brownovy zakřivené pyramidy. V~takovém případě totiž řešení musí splňovat $\sigma_1=\sigma_2=\sigma_3=\sigma_t$, odkud snadno odvodíme nutné a postačující podmínky pro tento typ návratu. Pokud návrat do vrcholu nenastává, testujeme, zda-li nenastane návrat na hladkou plochu tak, že najdeme pomocí Newtonovy metody kandidáta na neznámé  napětí $\mbf\sigma$. Pokud příslušná hlavní napětí splňují $\sigma_1\geq\sigma_2\geq\sigma_3$, potom je $\mbf\sigma$ skutečně hledaným řešením. Pokud tomu tak není, rozhodneme, zda-li nastane návrat na levou či pravou hranu podle toho, která z nerovností v~$\sigma_1\geq\sigma_2\geq\sigma_3$ neplatí. V~našem navazujícím výzkumu plánujeme ukázat, že tento ad-hoc algoritmus je dobře definovaný a spolehlivý, a to na základě subdiferenciálního přístupu.

Ke konstrukci zobecněné derivace konstitutivního operátoru $T$ lze použít tzv. vlastní projekce a jejich derivace, viz \cite{neto_computational_2008}. Pro stručnost budeme uvažovat pouze případ, kdy testovací tenzor malých deformací $\mbf\varepsilon^{tr}$ má tři různá vlastní čísla, tj. platí $\varepsilon_1^{tr}>\varepsilon_2^{tr}>\varepsilon_3^{tr}$. Potom spektrální rozklad $\mbf\varepsilon^{tr}$ splňuje
$$\mbf \varepsilon^{tr}=\sum_{i=1}^3\varepsilon^{tr}_i\mathbf e_i^{tr}\otimes\mathbf e_i^{tr},\quad \mathbf e_i^{tr}\otimes\mathbf e_i^{tr}=\mbf E_i^{tr}=\frac{(\mbf\varepsilon^{tr}-\varepsilon^{tr}_j\tn I)(\mbf\varepsilon^{tr}-\varepsilon^{tr}_k\tn I)}{(\varepsilon^{tr}_i-\varepsilon^{tr}_j)(\varepsilon^{tr}_i-\varepsilon^{tr}_k)},\quad i=1,2,3.$$
Říkáme, že tenzory $\mbf E_1^{tr}$, $\mbf E_2^{tr}$ a $\mbf E_3^{tr}$ jsou \textit{vlastními projekcemi} tenzoru $\mbf\varepsilon^{tr}$. Pokud uvažujeme vlastní čísla $\varepsilon_1^{tr}>\varepsilon_2^{tr}>\varepsilon_3^{tr}$ jako funkce závisející na tenzoru $\mbf\varepsilon^{tr}$, potom jejich derivace splňují $D\varepsilon^{tr}_i(\mbf\varepsilon^{tr})=\mbf E_i^{tr}$, $i=1,2,3$. Dále, derivace tenorových funkcí $\mbf E_i^{tr}$ jsou tenzory čtvrtého řádu a mají následující tvar:
\begin{eqnarray*}
	D\tn E_i^{tr}(\mbf{\varepsilon}^{tr})&=&\frac{\mathrm{D}((\mbf\varepsilon^{tr})^2)-(\varepsilon^{tr}_j+\varepsilon^{tr}_k)\tn{I} -(2\varepsilon^{tr}_i-\varepsilon^{tr}_j-\varepsilon^{tr}_k)\tn E_i^{tr}\otimes\tn E_i^{tr}}{(\varepsilon^{tr}_i-\varepsilon^{tr}_j)(\varepsilon^{tr}_i-\varepsilon^{tr}_k)}\\[2mm]
	&&-\frac{(\varepsilon^{tr}_j-\varepsilon^{tr}_k)[\tn E_j^{tr}\otimes\tn E_j^{tr}-\tn E_k^{tr}\otimes\tn E_k^{tr}]}{(\varepsilon^{tr}_i-\varepsilon^{tr}_j)(\varepsilon^{tr}_i-\varepsilon^{tr}_k)},
\end{eqnarray*}
kde $i\neq j\neq k\neq i$.

Položíme-li $\mbf\varepsilon^{tr}:=\mbf\varepsilon_k- \mbf \varepsilon_{k-1}^p+t_k\mbf \varepsilon_0/t_{max}$ a předpokládáme-li, že $\varepsilon_1^{tr}>\varepsilon_2^{tr}>\varepsilon_3^{tr}$, potom můžeme upřesnit definici konstitutivního operátoru $T$ a nejít jeho zobecněnou derivaci následovně, viz \cite{sysala_subdifferential-based_2017}:
$$\mbf\sigma_k=T\left(\tn{\varepsilon}^{tr}\right)=\sum_{i=1}^3\sigma_i(\mbf\varepsilon^{tr})\tn E_i^{tr},\quad \mathrm{D}  T\left(\mbf\varepsilon^{tr}\right)=\sum_{i=1}^3\left[\sigma_i(\mbf\varepsilon^{tr})D\tn E_i^{tr}+\mbf E_i^{tr}\otimes \mathrm{D}\sigma_i(\mbf\varepsilon^{tr})\right].$$
Zde $\sigma_1,\sigma_2,\sigma_3$ jsou vlastní čísla hledaného tenzoru napětí $\mbf\sigma_k$, která závisejí pouze na $\mbf\varepsilon^{tr}$, a $\mathrm{D}\sigma_i$ značí zobecněné derivace $\sigma_i$, $i=1,2,3$, které nalezneme podle známých vět o~derivacích implicitních funkcí.


\subsubsection{Numerický příklad}
\label{sec_example}

V této části porovnáme možnosti pružných a pružně-perfektně plastických modelů predikovat zóny poškození. Pro účely porovnání uvažujeme úlohu rovinné deformace inspirovanou případovou studii TSX experimentu vykonaného v~hloubce okolo 500 metrů v~URL, Kanada.

\begin{figure}
	\begin{center}	
		\includegraphics[width=0.4\textwidth]{graphics/mesh.pdf}	
		\includegraphics[width=0.45\textwidth]{graphics/stress_init.pdf}	
	\end{center}
	\caption{Geometrie a síť úlohy rovinné deformace. Rozměry jsou v metrech.}
	\label{fig_scheme_TSX}
\end{figure}

Geometrie úlohy a použitá konečně-prvková síť jsou znázorněny na \figref{fig_scheme_TSX}. Konkrétně, uvažujeme tunel s eliptickým profilem s průměry os 4.375 a 3.5 metrů. Počáteční tenzor napětí $\mbf\sigma_0$ je předpokládán konstantní v~celé oblasti a jeho nenulové komponenty mají následující velikosti: $\sigma_{0,1}=-45$ MPa, $\sigma_{0,2}=-11$ MPa a $\sigma_{0,3}=-60$ MPa. Největší hlavní napětí $\sigma_{0,3}$ ($\sigma_{0,3}$) působí ve směru osy tunelu a je zahrnuto do  modelu prostřednictvím Mohr-Coulombova plastického kritéria, ačkoliv je řešena 2D úloha. Zbývající hlavní napětí jsou znázorněna na \figref{fig_scheme_TSX}. Ražba tunelu trvala 17 dní, čemuž odpovídá hodnota parametru $t_{max}$. V~rámci výpočtu rozdělíme časový interval na 17 výpočetních kroků. Materiálové parametry volíme následovně: $E=60$ GPa, $\nu=0.2$, $c=17$ MPa a $\phi=\psi=26^\circ$. Poznamenejme, že hodnoty pevnostních parametrů $c$ a $\phi$ jsou nižší než v~\cite{Rutqvist2009Modeling}, abychom mohli snáze vizualizovat rozdíl mezi pružným a pružně-plastickým přístupem. Pro výpočet použijeme 6-bodové trojúhelníkové konečné prvky (P2) a 7-bodovou numerickou kvadraturu na jednotlivých prvcích. Pro implementaci úloh jsme použili a rozšířili vlastní kódy v~Matlabu. Některé z nich jsou volně stažitelné, viz \cite{cermak_efficient_2019}. Také je veřejně dostupná verze těchto kódů pro programovací jazyk Python, viz \cite{beseda_fem-elastoplasticity_2022}.

Porovnání zón poškození vypočtených pro pružné a pružně-plastické modely je znázorněno na obrázku \ref{fig_zones}. Zóny poškození pro pružný model jsou tvořeny konečnými prvky, na kterých není splněno Mohr-Coulombovo kritérium. Vidíme, že jejich tvar je zaoblený vůči stěně tunelu. V~případě pružně-plastického modelu zóny poškození reprezentují konečné prvky, na kterých je nenulová plastická deformace. Tyto zóny mají tvar V-zářezu, který je často pozorován v~in-situ experimentech. 
\begin{figure}
	\begin{center}	
		\includegraphics[width=0.42\textwidth]{graphics/plast_zone_elast.pdf}	
		\includegraphics[width=0.42\textwidth]{graphics/plast_zone.pdf}	
	\end{center}
	\caption{Porovnání zón poškození pro pružný model (vlevo) a pružně-plastický model (vpravo).}
	\label{fig_zones}
\end{figure}

Na \figref{fig_stress} je znázorněno porovnání horizontálních napětí pro oba přístupy. Vidíme, že elastický přístup povoluje vyšší koncentrace napětí (okolo 100 MPa) na stropě a počvě tunelu, zatímco v~případě pružně-plastického modelu jsou tyto koncentrace nižší, pouze okolo 50 MPa.
\begin{figure}
	\begin{center}	
		\includegraphics[width=0.45\textwidth]{graphics/stress_horizontal_elast.pdf}	
		\includegraphics[width=0.45\textwidth]{graphics/stress_horizontal_plast.pdf}	
	\end{center}
	\caption{Porovnání horizontálních napětí pro pružný model (vlevo) a pružně-plastický model (vpravo). Výsledky jsou v~MPa.}
	\label{fig_stress}
\end{figure}

Obrázek \ref{fig_deform} znázorňuje porovnání celkových posunů a současně 300 krát zvětšený deformovaný tvar tělesa. Pro lineární pružnost pozorujeme kontrakci horninového masivu na stropě a počvě tunelu, zatímco pro pružně-plastický model je pozorována dilatance. Pro lepší vizualizaci tohoto zásadního rozdílu byl vytvořen \figref{fig_displacement}, na kterém jsou porovnávány vertikální posuny v~závislosti na čase. Vidíme, že oba modely dávají shodné výsledky (kontrakci) do desátého časového kroku. Potom ale plastická odezva začíná být silně nelineární a vede k~objemové dilatanci.

\begin{figure}
	\begin{center}	
		\includegraphics[width=0.45\textwidth]{graphics/deform_shape_elast.pdf}	
		\includegraphics[width=0.45\textwidth]{graphics/deform_shape_plast.pdf}	
	\end{center}
	\caption{Porovnání celkových posunů a deformovaných tvarů pro pružný model (vlevo) a pružně-plastický model (vpravo). Výsledky jsou v~metrech.}
	\label{fig_deform}
\end{figure}
\begin{figure}
	\begin{center}	
		\includegraphics[width=0.8\textwidth]{graphics/vertical_displ.pdf}	
	\end{center}
	\caption{Vertikální posuny (v metrech) na stropě tunelu v~závislosti na čase.}
	\label{fig_displacement}
\end{figure}




\subsection{Pokročilé numerické metody řešení poroelastických problémů}
\label{sec:advaced_mech}
Tato sekce navazuje na kapitolu \ref{sec:poroelasticity} a uvedeme zde tři pokročilé metody řešení poroelastických problémů v~kontinuální oblasti i v~oblasti s puklinovou sítí, kterými jsme se zabývali v~rámci řešení projektu ENDORSE.

\subsubsection{Předpodmiňovače pro poroelastické problémy}

V sekci \ref{sec:biot_diskretizace} jsme představili metodu fixed-stressed splitting pro řešení poroelastických problémů. Nyní stručně představíme alternativní iterační řešič, který je založen na předpodmíněné Krylovovské metodu (například GMRES nebo MINRES). Omezíme se přitom pouze na lineární poroelastické problémy uvažované v~kontinuální oblasti, tedy na řešení systému lineárních rovnic popsaných rovnicí \eqref{eq:hm_alg_lin}. Splitting metody i předpodmíněné Krylovovské metody mají společnou myšlenkou v~tom, že řešení svázaného systému rovnic pro mechaniku a proudění převádějí na řešení dílčích nezávislých úloh mechaniky a proudění.
	
Předpodmínění pro lineární soustavu se čtvercovou maticí $\tn K$ a pravou stranou $\vc b$
	\begin{equation}
		\tn K \vc x = \vc b
  \label{eq:K_system}
	\end{equation}
	lze chápat jako nalezení vhodné čtvercové invertovatelné matice $\tn P$, pomocí které systém \eqref{eq:K_system} přepíšeme na tvar
 \begin{equation}
		\tn P^{-1}\tn K \vc x = \tn P^{-1}\vc b
  \label{eq:K_system2}
	\end{equation}
a předpokládáme, že matice $\tn P^{-1}\tn K$ má lepší vlastnosti pro numerické řešení než původní matice $\tn K$, konkrétně, že matice $\tn P^{-1}\tn K$ má lepší číslo podmíněnosti než $\tn K$.

V případě lineárního poroelastického problému formulovaného v~posunutích, \sy{pórových} tlacích a rychlostem  má matice $\tn K$ přirozenou blokovou strukturu, která odpovídá jednotlivým polím,
	\begin{equation}
		\tn K = \begin{bmatrix}
			\tn A & & \tn{B}_u^{T} \\
			& \tn M & \tn{B}_v^{T} \\
			\tn{B}_u & \tn{B}_v & -\tn D
		\end{bmatrix}, \quad \vc x = \begin{bmatrix}
		\vc u \\ \vc v \\ \vc p
	\end{bmatrix},
	\end{equation}
kde $\vc u$ je pole posunů, $\vc p$ \sy{pórových} tlaků a $\vc v$ je rychlostní pole. Diagonální členy odpovídají jednotlivým dílčím problémům, mimodiagonální vazbám mezi nimi. Matice $\tn K$ je navíc symetrická.

Vhodné předpodmínění úlohy lze provést pomocí blokově diagonální a symetrické matice
	\begin{equation}
		\tn P = \begin{bmatrix}
			\tn A + \gamma_1 \tn B_u^{T} \tn D^{-1} \tn B_u & & \\
			 & \gamma_2 \tn M + \gamma_3 \tn B_v^{T} \tn D^{-1} \tn B_v & \\
			 & & \gamma_4 \tn D
		\end{bmatrix}
	\end{equation}
kde $\gamma_1$, $\gamma_2$, $\gamma_3$ a $\gamma_4$ jsou kladné konstanty vhodně zvolené v~závislosti na materiálových parametrech a použité konstrukcí předpodmínění, viz např. \cite{luber_efficient_2022, hong_parameter-robust_2018}.

Blokově diagonální předpodmiňovač nám umožňuje v~rámci iterační metody MINRES řešit pouze menší systémy rovnic s maticemi $\tn A + \gamma_1 \tn B_u^{T} \tn D^{-1} \tn B_u$, $\gamma_2 M + \gamma_3 \tn B_v^{T} \tn D^{-1} \tn B_v$ a $\gamma_4 \tn D$, které reprezentují jednotlivá fyzikální pole. Výhodou tohoto přístupu je, že matici $\tn B_u^{T} \tn D^{-1} \tn B_u$ lze najít v~explicitním tvaru. Dále, blok $\tn A + \gamma_1 \tn B_u^{T} \tn D^{-1} \tn B_u$ lze řešit standardními metodami pro řešení elastických úloh, pro blok $\gamma_2 M + \gamma_3 \tn B_v^{T} \tn D^{-1} \tn B_v$ je možné použít standardní řešiče pro tzv. $H(\div)$ soustavy a blok $\gamma_4 \tn D$ lze řešit triviálně přímou inverzí. 


\subsubsection{Paralelní řešení kontaktních úloh metodou FETI}
\label{sec:feti}

Metoda rozložení oblasti FETI (Finite Element Tearing and Interconnecting) patří dnes k~nejefektivnějším metodám řešení obrovských soustav lineárních rovnic a úloh kvadratického programování. FETI rozdělí danou výpočetní oblast, respektive její síť, na $N_s$ nepřekrývajících se podoblastí $\Omega_s$ (obvykle tolik, kolik je k~dispozici výpočetních jader $N_c$). Používá se k~tomu například software ParMETIS, viz \cite{Karypis2000ParMetis}. Obrázek \ref{fig:feti} ilustruje rozložení oblasti pomocí FETI metody.
\begin{figure}
\begin{center}
\includegraphics[width=\textwidth]{graphics/tfeti2d.eps}
\caption{Ilustrace nepřekrývajícího rozložení čtvercové oblasti na čtyři podoblasti. Lagrangeovy multiplikátory $\lambda$ představují síly lepící podoblasti, respektive Dirichletovy okrajové podmínky.}
\label{fig:feti}
\end{center}
\end{figure}

Každá podoblast $\Omega_s$ je diskretizována metodou konečných prvků, čímž získáme matici tuhosti $\tn K^s$, její jádro $\tn R^s$ a zatěžující vektor $\vc f^s$ $s=1,\,\ldots,\,N_S$, ze kterých se sestaví globální distribuované objekty
%
\begin{flalign*}
\tn K &= \mbox{diag}(\tn K^1,\,\ldots,\,\tn K^{N_S}),\\
\tn R &= \mbox{diag}(\tn R^1,\,\ldots,\,\tn R^{N_S}),\\
\vc f &=[(\vc f^1)^T,\,\ldots,\,(\vc f^{N_S})^T]^T,
\end{flalign*}
kde diag značí blokově-diagonální matici obsahující na diagonále bloky uvedené v~závorkách.
Následně je třeba zajistit lepení těchto podoblastí (ztotožnění příslušných neznámých na rozhraních podoblastí) pomocí Booleanovské matice $\tn  B_E$ ($E=$ Equalities). Obdobně se dají zapsat podmínky nepronikání pomocí matice $\tn  B_I$ ($I=$ Inequalities). Dirichletovy okrajové podmínky se vynutí standardní modifikací matic $\tn  K^s$ a vektorů $\vc  f^s$ těch podoblastí, do kterých zasahuje tato podmínka. Výše popsaný přístup je tzv. FETI-1 varianta. 
Nevýhodou FETI-1 je, že jednotlivé matice $\tn  K^s$ mají různá jádra $\tn  R^s$, která mohou být obtížně spočitatelná. Proto preferujeme variantu metody Total-FETI (TFETI), kde předepisujeme Dirichletovy okrajové podmínek jako lineární rovnostní omezení pomocí matice $\tn  B_D$. Řádky matice $\tn B_D$ jsou typicky pouze přidány do matice $\tn B_E$. Výhodou je snadnější a stabilnější implementace, neboť všechny podoblasti mají stejná jádra $\tn  R^s$, která jsou numericky spočitatelná, viz \cite{gosselet_non-overlapping_2006} nebo sestavitelná z informace o~síti, viz \cite{dostal_total_2006}.
Matice $\tn  B_I$ je prázdná pro lineární úlohy (bez kontaktu). %Místo lepení případně nepronikání můžeme předepisovat požadovanou (minimální) vzdálenost těles skrze vektor $\vc {c}$.


Implementace FETI metod v~knihovně PERMON je následující: po získání globální primárních dat $\tn K$, $\vc f$,
$\tn B_E$, $\tn  B_{I}$ a $\tn R$, zavolá modul PermonQP specifickou sekvenci QP transformací za účelem získání QP úlohy s jednoduchým nerovnostním omezením a lineárním rovnostním omezením vyřešitelnou pomocí $QPSSolve()$ funkce.

Z matematického úhlu pohledu, volané QP transformace (QPT) implementují následující modifikace. Původní úloha 
\begin{equation}
 \argmin_{{\vc  u}}\frac{1}{2}\vc  u^{T}\tn K \vc u-\vc  f^{T}\vc  u\;\;\text{s.t.}\;\;\tn  B_{I}\vc  u\leq0\;\;\text{a}\;\;\tn  B_{E}\vc  u=0,\label{ekv_16-1}
\end{equation}
je transformována na duální pomocí $QPTDualize()$
\begin{equation}
 \argmin_{{\vc  \lambda}}\frac{1}{2}\vc  \lambda^{T}\tn F \vc \lambda-\vc  \lambda^{T}\vc  d\;\;\text{s.t.}\;\;\vc  \lambda_{I}\geq0\:\;\text{a}\;\:\tn G \vc \lambda=e,\label{ekv_24}
\end{equation}
kde vektor $\vc \lambda$ reprezentuje Lagrangeovy multiplikátory ilustrované na \figref{fig:feti},
\begin{equation}
\tn  F=\tn  B\tn  K^{\dagger}\tn  B^{T},\,\tn  G=\tn  R^{T}\tn  B^{T},\,\vc  d=\tn  B\tn  K^{\dagger}\vc  f,\,\vc{e}=\tn  R^{T}\vc  f \label{eq:FETIsymbols},
\end{equation}
přičemž $\tn  K^{\dagger}$ značí zobecněnou inverzi matice $\tn  K$ splňující $\tn  K\tn  K^{\dagger}\tn  K=\tn  K$.
Matice s omezeními 
$\tn  B=\left[\begin{smallmatrix}\tn  B_{I}^T & \tn  B_{E}^T\end{smallmatrix}\right]^T$
může být sestavena tak, že bude mít plnou hodnost a  Hessián $\tn F$ je pak pozitivně definitní matice s příznivě rozloženým spektrem pro metodu sdružených gradientů (CG). Z pohledu paralelizace je důležité, že akce $\tn  K^{\dagger}$ je pro každou podoblast lokální.

Řešení $\vc  u$ může být vyčísleno podle formule
\begin{equation}
\vc  u=\tn  K^{\dagger}(\vc  f-\tn  B^{T}\vc  \lambda)+\tn  R\vc{\alpha},\label{eq3}
\end{equation}
kde
\[
\vc{\alpha}=-(\tn  R^{T}\widetilde{\tn  B}^{T}\widetilde{\tn  B}\tn  R)^{-1}\tn  R^{T}\widetilde{\tn  B}^{T}\widetilde{\tn  B}\tn  K^{\dagger}(\vc  f-\tn  B^{T}\vc  \lambda)
\]
značí vektor amplitud tuhých módů určující příspěvek jádra $\tn  R\vc{\alpha}$ do $\vc  u$. Matice $\tilde{\tn  B}$ je definována jako $\tilde{\tn  B}=\left[\begin{smallmatrix}\tilde{\tn  B_{I}}^T & \tn  B_{E}^T\end{smallmatrix}\right]^T$ s $\tilde{\tn  B_{I}}$ tvořenou řádky $\tn  B_{I}$, které odpovídají aktivním omezením.

Minimalizační úloha na podmnožině afinního prostoru je transformována na úlohu na vektorovém podprostoru pomocí $\widetilde{\vc  \lambda}$ z oboru hodnot $\tn  G^T$ splňujícímu $\tn  G\widetilde{\vc  \lambda}=\vc{e}$, přičemž řešení má tvar $\vc  \lambda=\widehat{\vc  \lambda}+\widetilde{\vc  \lambda}$, kde $\widehat{\vc  \lambda}$ je z jádra matice $\tn  G$.

Úloha, kterou obdržíme pomocí $QPTHomogenizeEq()$ vypadá následovně
\begin{equation}
 \argmin_{\widehat{\vc  \lambda}}\frac{1}{2}\widehat{\vc  \lambda}^{T}\tn  F\widehat{\vc  \lambda}-\widehat{\vc  \lambda}^{T}(\vc  d-\tn  F\widetilde{\vc  \lambda})\;\;\text{s.t.}\;\;\widehat{\vc  \lambda}_{I}\geq-\widetilde{\vc  \lambda_{I}}\:\;\text{a}\;\:\tn  G\widehat{\vc  \lambda}=0.\label{ekv_25}
\end{equation}

Další vylepšení je založeno na vynucení homogenizovaného rovnostního omezení projektorem $\tn{P=I-Q}$ na jádro $\tn  G$, kde $\tn {Q}=\tn  G^{T}(\tn{GG}^{T})^{-1}\tn  G$ je ortogonální projektor na obor hodnot $\tn  G^{T}$: $\text{Im}\tn{P}=\text{Ker}\tn  G,\text{\,Im}\tn{Q}=\text{Im}\tn  G^{T}$.


Vyčíslení $(\tn{GG}^{T})^{-1}$, tj. řešení soustavy lineárních rovnic
\begin{equation} \label{ekvCP}
\tn{GG}^T\vc{x}=\vc{y},
\end{equation}
se nazývá hrubý problém.
Takto modifikovaná úloha získaná z \eqref{ekv_25} pomocí transformace $QPTEnforceEqByProjector()$ má následující tvar:
\begin{equation}
 \argmin_{\widehat{\vc  \lambda}}\frac{1}{2}\widehat{\vc  \lambda}^{T}\tn{PFP}\widehat{\vc  \lambda}-\widehat{\vc  \lambda}^{T}\tn{P} \vc d\;\;\text{s.t.}\;\;\widehat{\vc  \lambda}_{I}\geq-\widetilde{\vc  \lambda_{I}}\:\;\text{and}\;\:\tn  G\widehat{\vc  \lambda}=0.\label{ekv_26}
\end{equation}
%
Více detailů lze najít např. v~\cite{kozubek_solving_2016}.



\subsubsection{Nelineární efekty na puklinách}
\label{sec:nonlinear_fr_effects}

Při modelování hydro-mechanických procesů v puklinách hrají významnou roli nelineární jevy.
Deformace matrice, která se v měřítku lokality jeví malá, a tedy popsatelná pomocí lineární pružnosti, může způsobovat velkou relativní změnu rozevření puklin, k~jejímuž popisu je nutná nelineární mechanická teorie.
Tento fakt je částečně zachycen díky kontaktním podmínkám popsaným v~kapitole \ref{sec:biot_dfm}.
Zde uvedeme další jevy, jež mohou přispět k~realističtějšímu modelování hydromechaniky v oblastech s puklinami.

Změna objemu pukliny vyvolaná deformací matrice má vliv na hydraulickou vodivost.
Běžně používaným vztahem je tzv. kubický zákon (viz \cite{Snow1969Anisotropic})
\eq{ k_f(\uu) = \frac{\eta \varrho g}{12\mu} a(\uu)^2, }
kde $k_f\;\bunit{m\cdot s^{-1}}$ je hydraulická vodivost pukliny, $\eta\in(0,1]$ $[-]$ je koeficient drsnosti, $g\;\bunit{m\cdot s^{-2}}$ je gravitační zrychlení, $a(\uu) = \delta + \uu_m^+\cdot\vnu^+ + \uu_m^-\cdot\vnu^-\;\bunit{m}$ je šířka (rozevření) deformované pukliny, $\varrho\;\bunit{kg\cdot m^{-3}}$ a $\mu\;\bunit{Pa\cdot s}$ je hustota, resp. viskozita kapaliny.

Drsnost povrchu pukliny způsobuje další nelineární jev zvaný smyková dilatace.
Jejím principem je, že vlivem nerovností vede mechanické namáhání v~tečném směru k~rozevírání pukliny ve směru příčném. Matematicky lze smykovou dilataci chápat jako změnu kontaktní délky $\delta_{min}$ $\bunit{m}$ (minimální šířky pukliny, při níž dochází ke kontaktu) v~závislosti na relativním tečném posunu.
Kontaktní podmínka $\eqref{eq:contact}_2$ pak má tvar
\eq{ a(\uu) \ge \delta_{min}(\uu), }
kde $\delta_{min}$ je dáno profilem drsnosti, např.
\eq{ \delta_{min}(\uu) = \min\{|\uu_t|\operatorname{tg}\varphi, h\}. }
Zde $\uu_t = (\uu_m^+-\uu_m^-)_t$ $\bunit{m}$ je relativní tečné posunutí protilehlých stěn pukliny, $\varphi$ $\bunit{rad}$ a $h$ $\bunit m$ je úhel, resp. výška drsnosti (viz Obr. \ref{fig:rough_frac}).
%\sy{[Pozn. od EH: Doporučuji detailněji vysvětlit, o jaký úhel se jedná- třeba s pomocí obrázku.]}


\begin{figure}
\centering
\includegraphics[width=0.5\textwidth]{graphics/rough_frac.pdf}
\caption{Profil drsnosti pukliny a jeho parametrizace v modelu smykové dilatace.}
\label{fig:rough_frac}
\end{figure}


Obě zmíněné nelinearity, vodivost závislou na rozevření a smykovou dilataci, lze začlenit do numerického schématu popsaného v~sekci \ref{sec:biot_diskretizace}, které je založeno na iteračním štěpení hydrauliky a mechaniky.

\subsection{Homogenizace hydraulické vodivosti}
\label{sec:hm_homogenizace}

Vliv EDZ na hydraulickou vodivost lze nejjednodušším způsobem popsat pomocí empirických vztahů (\cite{Rutqvist2009Modeling,Snow1969Anisotropic}). Například v kapitole \ref{sec:TSX} jsme uvedli vztah \eqref{eqn:nonlinear_cond_exp}.
Tento postup ale nezohledňuje heterogenitu a anizotropii prostředí.
Experimentální i numerické studie (např. \cite{Min2004Stressdependent,Bidgoli2013Numerical}) však potvrzují, že tenzor vodivosti je v porušené zóně značně anizotropní, přičemž zásadní roli hraje (náhodná) mikrostruktura poruch.
V této kapitole popíšeme přístup pro určení vlivu puklinové sítě a mechanického namáhání na hydraulickou vodivost, založený na numerické homogenizaci.
Použitím tohoto postupu je možné dále zpřesnit odhad hydraulických parametrů EDZ.

Uvažujeme jednoduchou geometrickou oblast (čtverec, krychle) reprezentující vzorek porušené horniny, která obsahuje síť náhodných diskrétních puklin. Pro názornost se zde omezíme na rovinnou oblast $\Omega=(0,1)\times(0,1)$, tedy čtverec o~straně 1 metr.
Na této oblasti řešíme $n$ hydraulických úloh (případně poroelastických úloh, chceme-li zahrnout vliv mechaniky na vodivost). Jako okrajovou podmínku pro hydraulický model volíme pro $i=1,...,n$ předepsaný tlak: 

% \sy{[Standa: Nechybí zde nějaká suma přes index $i$, případně neměl by index i být i na levé straně? Takhle mi to moc nedává smysl.]}
% \js{[JS: Pokusil jsem se naformulovat srozumitelněji.]}

$$ p_i(\vc x) = \eps \vc x\cdot\vc t_i, \quad \mbox{kde } \vc t_i=\begin{bmatrix}\cos\omega_i\\\sin\omega_i\end{bmatrix}, ~\omega_i=i\tfrac{2\pi}{n}. $$
Vektor $\eps\vc t_i$ představuje průměrný tlakový gradient v oblasti.
Okrajové podmínky pro mechaniku stanovíme dle předpokládaného mechanického namáhání a časový interval volíme tak, aby na jeho konci došlo k ustálení HM procesů.
Pro $i=1,...,n$ pak definujeme průměrné pole rychlosti 
$$ \eps\vc v_i = \frac1{|\Omega|}\sum_{E\in\mathcal T}|E|\vc v_E, $$
kde $\vc v_E$ je rychlost na elementu $E$ v koncovém čase.
Tyto hodnoty použijeme pro fitování homogenizovaného tenzoru hydraulické vodivosti
$$ \tilde{\tn K}=\begin{bmatrix}k_{11}&k_{12}\\k_{12}&k_{22}\end{bmatrix} $$
metodou nejmenších čtverců na Darcyův zákon
$$ -\tilde{\tn K}\vc t_i = \vc v_i,~i=1,...,n. $$
Takto získaný tenzor $\tilde{\tn K}$ pak již zohledňuje vodivost intaktní horniny, vodivost puklin i mechanické vlivy.

Na \figref{fig:hm_homogenizace} je příklad numerické homogenizace za použití poroelastické úlohy s vodivostí puklin závislou na rozevření.
Pro režim bez zatížení a pro smykové, resp. tlakové zatížení je zobrazen rozdíl vodivosti puklin oproti klidovému stavu s nulovým tlakovým gradientem.
Dále jsou vyobrazena vlastní čísla a vlastní vektory homogenizovaného tenzoru hydraulické vodivosti. Je vidět silný vliv zatížení na anizotropii $\tilde{\tn K}$.
\begin{figure}
\centering
\includegraphics[width=\textwidth]{graphics/hm-homogenizace.png}
\caption{Homogenizace tenzoru hydraulické vodivosti v oblasti s náhodnými puklinami. Elipsy znázorňují vlastní čísla a vlastní vektory tenzoru $\tilde{\tn K}$ pro různé případy smykového, resp. tlakového zatížení.}
\label{fig:hm_homogenizace}
\end{figure}

% \sy{[Standa: Numerický příklad, jeho vstupy a výstupy bych doporučil popsat podrobněji.]}
% \js{[JS: Stačí takto nebo ještě podrobněji?]}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
\section{Návrh in-situ experimentu za účelem ověření navržené metodiky}
\label{sec:Bukov}

%\sy{[Standa: V této kapitole jsem použil loňský text s návrhem experimentu a prvotně jej napojil na metodiku, aby byl text aspoň trochu kompaktní. Zkrátil jsem úvod a vyhodil část věnovanou experimentu TSX, protože je již popsaný v~šesté kapitole. Naopak do šesté kapitoly jsem přidal více detailů o~TSX. Dále, v~této kapitole jsem přepsal EIZ na EDZ, neboť v~rámci metodiky jsme tyto pojmy sloučili. Návrh experimentu doporučuji psát více obecně, protože jeho vykonání na PVP Bukov se jeví jako problémové kvůli přílišnému odvodnění. Prosím o~revizi této kapitoly, protože mnoha věcem příliš nerozumím.]}

Navrženou metodiku lze pouze částečně validovat vůči existujícím experimentům, především vůči dílčím měřením experimentu TSX (URL Kanada) a měřením v rámci experimentu FEBEX (Grimsel, Švýcarsko). Proto se v~této kapitole zaměříme na návrh experimentu, který by umožnil tuto metodiku otestovat a případně také certifikovat. Hlavním principem navrhovaného experimentu je kontinuální měření pórového tlaku v~různé vzdálenosti od stěny díla v~průběhu ražby a následně po dobu relaxace horninového masivu. Pomocí inverzních metod pak bude určena permeabilita v~EDZ včetně anizotropie, heterogenity a parametrů modelu závislosti permeability na změnách tenzoru napětí v hornině. Alternativně bude permeabilita predikována též pomocí cenově dostupnějších, ale méně relevantních geofyzikálních měření. Experiment by mohl být uskutečněn v~rámci PVP Bukov či jiné podzemní laboratoře.

Experiment TSX byl již popsán v hlavní části metodiky, pro úplnost uvádíme v části \ref{subsec:FEBEX} stručný popis měření tlaků v rámci experimentu FEBEX. Následně definujeme cíle navrhovaného experimentu v~části \ref{sec:cile} a dále v~části \ref{sec:navrh} popíšeme postup prací při realizaci experimentu. V~části \ref{sec:predikce} popíšeme návrh modelářských prací, které by měly předcházet finálnímu návrh experimentu, a orientačně porovnáme klíčové parametry hornin na PVP Bukov vůči podobným experimentům v~rámci projektů TSX (URL) a FEBEX (Grimsel).  V~části \ref{sec:analyza} pak navrhneme postup při zpracování dat z experimentu a validaci metodiky predikce indikátorů EIZ. V~poslední části \ref{sec:swot} předkládáme stručnou SWOT analýzu návrhu.

\subsection{Experiment FEBEX}
\label{subsec:FEBEX}

Experiment FEBEX (1994--2000), provedený na lokalitě Grimsel ve Švýcarsku, byl zaměřen na dlouhodobé studium teplotních hydraulických mechanických a transportních (THMC) procesů při interakci bentonitové bariéry s okolní horninou v~konceptu horizontálního ukládání. V~rámci rozsáhlého experimentu bylo též realizováno kontinuální měření \sy{pórového} tlaků při ražbě ukládacího tunelu o~průměru 2.27~m a délce 70~m. Před vlastní ražbou byl realizován paralelní vrt ve vzdálenosti asi 2--3~m od stěny ukládacího vrtu. V~paralelním vrtu byly zapakrovány úseky P3 (62--74~m) a P4 (přibližně 48--62~m). V~těchto úsecích byl měřen tlak během následující ražby. Ražba byla realizována štítovým razícím strojem (TBM) s minimálním poškozením vlivem ražby (odhadováno do vzdálenosti 2~cm od stěny díla). Průběh ražby, umístění měřených intervalů a pozorovaný průběh tlaků shrnuje  
\figref{fig:FEBEX}. 

\begin{figure}[h!]
\includegraphics[width=\textwidth]{graphics/FEBEX_drilling_pressure_outlined.pdf}
\caption{Průběh měřených tlaků při ražbě tunelu FEBEX., zdroj \cite{rutqvist_analyses_2004}.}
\label{fig:FEBEX}
\end{figure}

Měření zřetelně odráží vliv ražby, avšak následný pokles \sy{pórových} tlaků je výrazně rychlejší (v řádu desítek hodin) nežli v~případě experimentu TSX. To je zřejmě způsobeno řádově větší permeabilitou (cca $10^{-10}\,\unit{m\cdot s^{-1}}$ pro okolí intervalu P4). Překvapující je fakt, že při měření na boku tunelu dochází k~nárůstu \sy{pórového} tlaku, ačkoliv by to mělo naopak docházet k~jeho poklesu.
Výsledky inverzního modelování s využitím 2D i 3D modelů vedou k~závěru, že lze dosáhnout dobré shody modelu s měřením za předpokladu otočení tenzoru počátečního napětí o~přibližně $40\degree$.  Vzhledem k~velmi komplikované geologické situaci nelze tuto možnost vyloučit. Měřený tenzor počátečního napětí má hlavní komponenty 30, 15, 10~MPa  s maximálním napětím pod úhlem $45\degree$ od směru vrtu a minimálním napětím ve svislém směru, nejlepší shody s měřením bylo dosaženo pro hlavní komponenty 30, 15, 22~MPa, avšak s maximálním napětím ve směru vrtu a zvýšenou hodnotou 22~MPa  pro vertikální napětí. 

\subsection{Cíle experimentu}
\label{sec:cile}
Z hlediska bezpečnosti je klíčovým parametrem EDZ tenzor permeability $\perm$, respektive tenzor hydraulické vodivosti
\[
    \tn K = \perm \frac{\rho g}{\eta},
\]
kde $\rho$, $\eta$ a $g$ značí postupně hustotu, viskozitu vody a gravitační zrychlení. Tenzor permeability $\perm$ je třeba predikovat modelem, jelikož je závislý na tenzoru deformace způsobené změnami napětí, ke kterým bude docházet i po uzavření úložiště v~důsledku opětovného zatopení a v~menší míře i v~důsledku bobtnacího tlaku bentonitu. Pro posouzení bezpečnosti jednotlivých úložných míst z hlediska EDZ je třeba určit i variabilitu permeability v~prostoru, zejména ve směru úložné chodby. Předpokládáme, že model změny permeability v~důsledku deformace, je nezávislý na poloze, a proto uvažujeme permeabilitu v~separovaném tvaru: 
\[
    \perm(\vc x, \strain) = \perm'(\strain)\perm_0(l),
\]
kde $\perm_0(l)$ popisuje počáteční heterogenitu horniny pro staničení $l$ a $\perm'$ je model závislosti změny permeability na tenzoru deformace $\strain$ (případně tenzoru napětí). Deformace se pak mění zejména v~rovině kolmé na osu díla, změnu ve směru chodby by bylo vhodné uvažovat při výrazné variabilitě mechanických parametrů horniny ve směru chodby. 

Hlavní cíle experimentu jsou:
\begin{itemize}
    \item validace různých modelů definujících funkci $\perm'(\strain)$,
    \item validace postupů pro kalibrace parametrů těchto modelů,
    \item validaci postupů pro kalibraci prostorové variability permeability: $\perm_0(l)$,
    \item predikce transportních parametrů EDZ a vliv uvedených změn permeability na navržený indikátor bezpečnosti.
\end{itemize}
Pro každý z těchto čtyř cílů může být uvažováno více variant (modelů resp. postupů), součástí validace pak bude i výběr vhodné varianty.  
Možné vedlejší cíle zahrnují:
\begin{itemize}
     \item testování kvality predikce pouze na základě dat z průzkumného vrtu (před ražbou),
    \item ověření korelace mezi permeabilitou a výsledky geofyzikálních metod, v~geologických podmínkách PVP Bukov.
    \item ověření predikce indikátorů bezpečnosti EDZ na základě geofyzikálních měření
    \item porovnání odezvy pórového tlaku na ražbu v~geologických podmínkách PVP Bukov vůči uvedeným zahraničním experimentům.
\end{itemize}
Vedlejší cíle nejsou chápány jako méně důležité, ale spíše zatížené větší výzkumnou nejistotou. Jejich nesplnění 
by proto nemělo být chápáno jako neúspěch experimentu.

Po uzavření úložiště v~delším časovém horizontu je EDZ vystavena i dalším tepelným a chemickým vlivům, které mohou ovlivňovat její permeabilitu. Tyto vlivy pro jednoduchost nebudeme uvažovat v~plánovaném experimentu. Studované metody však umožní identifikaci klíčových parametrů hydro-mechanických procesů, které lze použít i v~případě zahrnutí teploty. Doplňkové parametry pak jsou zejména: teplotní roztažnost, teplotní vodivost a tepelná kapacita horniny, které lze poměrně dobře měřit laboratorně. 

\subsection{Koncept experimentu}
\label{sec:navrh}
V této části podrobněji představíme navrhovaný experiment včetně  doprovodných měření. Popíšeme navrhovanou posloupnost činností při realizaci experimentu a orientační množství a uspořádání realizovaných měření. Dále navrhneme přidružená měření, která by bylo vhodné realizovat pro lepší porozumění tvorbě EDZ a podpoření vedlejších cílů, ale nejsou nezbytná pro dosažení hlavních cílů experimentu. Nakonec shrneme klíčové požadavky na realizaci prací pro dosažení cílů experimentu. 

\subsubsection{Průběh experimentu}
Schematický návrh realizace experimentu je následující:
\begin{enumerate}
%    \sy{volba lokality...}
% JB: není součást experimentu, ale je popsáno v části 4.3, upřesnil jsem to tam.
    \item Před ražbou experimentálního tunelu budou v~nadloží a na boku tunelu vyvrtány dva měřicí vrty ve vzdálenosti cca 1m od plánované stěny díla. Pro dosažení lepších výsledků následných inverzních úloh by bylo vhodné přidat další dvojici vrtů ve vzdálenosti cca 2.5~m od stěny díla.
    \item Na vybraných vrtech bude provedeno optické mapování pro získání puklinových statistik. Ty budou použity pro tvorbu mikroškálového modelu.
    \item Bude provedena minimálně jedna tlaková zkouška pro stanovení permeability neporušené horniny. Vzhledem k možným lokálním vlivům doporučujeme minimálně tři měření na různých místech se stejnou geologickou stavbou.
    \item Bude provedeno měření tenzoru počátečního napětí více různými metodami. Vhodné je zejména stanovení počátečního napětí pomocí konvergenčních měření případně měření metodou hydrofrac neboť identifikují efektivní počáteční napětí na škálác srovnatelných s rozměry tunelu a EDZ. 
    %\jb{Je reálné?}
    \item Měřící vrty budou osazeny tlakovými senzory v~zapakrovaných úsecích (délka cca 1~m) s roztečí  2--6~m s kabeláží pro kontinuální měření.  
    \item Bude provedeno kontrolní měření pórových tlaků v~klidovém stavu po dobu min. 1 měsíce.
    \item Bude realizována ražba experimentální chodby v~návrhové délce odpovídající alespoň třem roztečím mezi úložnými pozicemi. %\jb{předpokládanou rozteč kontejnerů se mi nepodařilo nalézt}. 
    Předpokládáme průměr chodby 3-5m, pro minimalizaci EDZ je vhodné volit kruhový resp. eliptický průřez. Pro přesnější analýzu bude vhodné znát přesně výsledný tvar chodby, např. pomocí laser skenu. Nutné je použití přesné a šetrné ražby pro minimalizaci EDZ, cílem je napodobit ražbu reálné úložné chodby.
    \item Kontinuální měření pórových tlaků bude probíhat po dobu minimálně jednoho roku od ražby. 
    \item Po ražbě bude na stěně experimentální chodby provedeno geologické mapování puklin. Bude realizováno geofyzikální měření na stěnách chodby metodami odporové tomografie, aktivní seismiky případně též měření georadarem. Pro tato měření je třeba mít k dispozici přesné prostorvé lokace pro řešení příslušných inverzních úloh.
    \item S odstupem cca 6 měsíců po ražbě budou realizovány 4 krátké vrty (cca 5~m) kolmo na boční stěnu díla a bude provedeno profilové měření permeability pomocí metody SEPPI a metody MVP.
    Měření bude opakováno cca 12 měsíců od ražby.
\end{enumerate}

Vlastní realizaci experimentu bude předcházet přípravná fáze zahrnující sběr dostupných dat a prediktivní modelování (viz kapitola \ref{sec:predikce}). Na základě této předběžné analýzy budou upřesněny parametry experimentu, zejména lokace, směr ražby, profil tunelu a pozice měření. 

\subsubsection{Vhodná doplňující měření a experimenty}
Pro zpřesnění interpretace hlavních měření a omezení rizika nedostatečné přesnosti navrhujeme realizovat též následující činnosti:
\begin{itemize}
    \item Konvergenční měření posunů na stěnách tunelu po jednotlivých odstřelech ražby za účelem určení počátečního napětí a jeho makroskopické variability pomocí inverzní analýzy. Konvergenční měření být v budoucnu nahrazeno pokročilou analýzou dat z laserového skenování (přesné měření délek mezi několika danými body je nahrazeno průměrováním přes velké množství méně přesných dat). Pro srovnání může být v hodné provádět ve stejných časových okamžicích oba typy měření. 
    \item Kontinuální měření změn napětí v~několika bodech podél experimentální chodby.  Z technických důvodů patrně v~samostatném paralelním vrtu.
    \item Monitoring akustických emisí pro upřesnění rozsahu poškození horniny.
    \item Ražba chodby rozdělená do dvou časově oddělených fází. Časový odstup je nutné stanovit na základě odhadů permeability a odpovídajících časů relaxace horniny. Vícefázová ražba by umožnila lépe studovat vliv ražby na vývoj pórového tlaku.
    \item Jádrové vrtání, stanovení \uv{Geological Strength Index} a triaxiální test pro určení parametrů plastického modelu: koheze, úhel vnitřního tření, maximální axiální tlak.
\end{itemize}
 
\subsubsection{Požadavky na správný průběh experimentu}
V této části shrneme okolnosti podstatné pro dosažení cílů experimentu.
\begin{itemize}
    \item Finálnímu návrhu experimentu by mělo předcházet jeho prediktivní modelování, viz následující kapitola \ref{sec:predikce},  pro posouzení citlivosti měření a optimalizaci volitelných parametrů uspořádání experimentu. Prediktivní modelování umožní: upřesnění lokace, orientace, profilu a rozměrů experimentální chodby, upřesnění polohy měřících vrtů a jejich případného úklonu, rozteč mezi tlakovými sondami, upřesnění postupu mapování trhlin.
    \item Dosavadní měření počátečního napětí na různých lokacích PVP Bukov ukazují na jeho značnou prostorovou heterogenitu. Podle expertního odhadu geologie na okolních chodbách a případných předchozích měření počátečních napětí bude vhodné zvolit dostatečné prostorové rozlišení pro měření počátečních napětí. Analýza časové posloupnosti laser skenů může umožnit volit plynule mezi přesností získaného pole počátečních napětí a jeho prostorovým rozlišením bez nutnosti optimalizace uspořádání měření. 
    \item Přesnou technickou realizaci měření tlaku je třeba konzultovat s odborníky s praktickou zkušeností.
    \item Pro všechny práce je nutné zajistit vzájemnou a otevřenou komunikaci mezi skupinami realizujícími plánovaní experiment, vlastní měření a interpretaci získaných měření.
    \item Minimálně po dobu šesti měsíců od ražby a ideálně až do ustálení měřených tlaků nesmí v~okolí cca 30~m probíhat žádná další ražba nebo tlaková zkouška, nebo jiné zásahy, které by mohly ovlivnit pórový tlak v~okolí tunelu.
    \item Měření SEPPI a MVP musí být realizováno tak, aby neovlivnilo kontinuální měření pórových tlaků.
\end{itemize}

\subsection{Prediktivní modelování}
\label{sec:predikce}

Před vlastní ražbou je vhodné provést prediktivní modelování experimentu, abychom na jedné straně zachytili požadované efekty (např. pozorovatelná tvorba EDZ), ale současně minimalizovali porušení horniny jako v~případě reálného úložiště. Ve spolupráci se SÚRAO bude nejprve stanoveno konkrétní umístění experimentální chodby, případně ve dvou variantách s odlišnou orientací vůči hlavním osám počátečního napětí. Následně budou shrnuta data o~horninách z již dostupných blízkých měření a na jejich základě budou pro vybrané varianty experimentu odhadnuty parametry horniny v~prostorách experimentu včetně jejich nejistoty. 

Na základě odhadnutých dat bude sestaven 2D model (řez kolmo na osu experimentální chodby) pro popis hydro-mechanických procesů v~důsledku ražby a predikci měřených pórových tlaků. 
Pro několik variant parametrů horniny bude provedena dopředná simulace měření tlaků, následně bude přičten časově korelovaný šum a bude provedena Bayesovská inverze. Rozptyl získaných aposteriorních rozdělení parametrů modelu odpovídá citlivosti měření na tyto parametry. Z možných variant pak bude vybrána ta, která zajistí dostatečnou citlivost na permeabilitu a její změny.

Pro vybrané varianty bude následně proveden výpočet pomocí 3D modelu se zahrnutím vlivu předepsané prostorové změny permeability, změny počátečního napětí ve směru experimentální chodby a náhodných prostorových fluktuací odpovídajících nejistotě dostupných dat. Pro simulované průběhy měřených tlaků bude opět provedena Bayesovská inverze (zde je šum již zahrnut v~prostorových fluktuacích) a získaná prostorově proměnná permeabilita bude porovnána vůči předepsané permeabilitě a vůči směrodatné odchylce aposteriorních rozdělení. Takto lze predikovat citlivost experimentu na permeabilitu se zahrnutím nežádoucího vlivu variability počátečního napětí. Pro stanovení počtu tlakových senzorů a tím prostorového rozlišení experimentu budou realizovány modelové inverzní úlohy s prostorovou změnou permeability na různých délkových škálách. Podobným způsobem bude pomocí 3D modelu testována a optimalizována citlivost parametrů závislosti permeability na deformaci. Zejména bude otestováno, zda uklonění vrtů vůči ose experimentální chodby může zvýšit tuto citlivost.
%JB: Zde si lze představit mnoho různých modelových výpočtů a část z nich bychom měli zkusit i
% v rámci ENDORSE.

\subsubsection{Orientační vlastnosti hornin na PVP Bukov a okolí}
V této části shrneme základní parametry hornin v~okolí PVP Bukov podle  \cite{Soucek2017Geotechnical}. Hydraulická vodivost se řádově pohybuje v~intervalu $10^{-13}$--$10^{-11}\, \unit{m\cdot s^{-1}}$, tedy je nižší než v~experimentu FEBEX a zároveň vyšší než v~experimentu TSX. Lze tedy očekávat \sy{ustálení pórových} tlaků po ražbě v~horizontu jednotek měsíců. Počáteční napětí bylo měřeno pomocí overcoringu (CCBO), hydrofrackingu (HF) a inverzní analýzou konvergenčních měření (INVGEM) s relativně velkými odchylkami způsobenými pravděpodobně značnou heterogenitou geologických podmínek. 
Vertikální hlavní napětí je uváděno v~rozsahu 20--30~MPa, horizontální hlavní napětí jsou přibližně v~intervalu 11--20~MPa s poměrem $1$ až $2$. Podle pozdějších měření metodou HF je rozsah horizontálních napětí 17--38~MPa s poměrem hlavních komponent 2--3, \cite{Bukovska2019Characterization}. Absolutní velikosti hlavních napětí jsou srovnatelné s experimentem FEBEX. Vzhledem k~velkým rozptylům nelze bez dalších dat a konzultací dělat žádné závěry o~poměrech hlavních napětí.

\subsection{Validační fáze}
\label{sec:analyza}

Po realizaci experimentu bude následovat validační fáze, ve které se na základě numerických simulací a inverzní analýzy naleznou vhodné materiálové parametry použitých hydro-mechanických modelů popisujících vznik EDZ a změnu hydraulické vodivosti. Inverzní úlohy budou realizovány pro čtyři sady dat:
\begin{itemize}
    \item data získaná před vlastní ražbou,
    \item data získaná před ražbou plus data z geofyzikálních měření a SEPPI,
    \item data získaná před ražbou plus data z kontinuálního měření tlaků,
    \item použití všech dostupných dat jako referenční popis EDZ.
\end{itemize}

První varianta napodobuje predikci před vlastní ražbou konkrétní úložné chodby. Druhá varianta modeluje teoreticky méně náročnou metodu charakterizace EDZ, bez paralelních vrtů, kdy je změna permeability stanovena pomocí lokálního profilování metodou SEPPI a prostorová variabilita odhadnuta pomocí geofyzikálních měření. Třetí varianta by měla poskytnout přesnější informaci o~změně permeability v~EDZ. 
Čtvrtá varianta bude použita pro posouzení věrohodnosti odhadů předchozích variant. Upřesnění validační fáze lze provést po zpřesnění návrhu celého experimentu.

\subsection{SWOT analýza}
\label{sec:swot}

\subsubsection{Silné stránky návrhu}
\begin{itemize}
    \item Návaznost na podobné experimenty v~zahraničí. Možnost uplatnění navržené metodiky a softwaru vyvinutého v~rámci TAČR projektu ENDORSE.
    \item Propojení in-situ experimentu s prediktivním modelováním, validačními numerickými simulacemi a Bayesovskou inverzní analýzou nejistot je velice komplexní přístup, který je přinejmenším v~českých podmínkách unikátní.
    \item Kontinuální měření signifikantních změn pórového tlaku umožňuje získat velké množství dat pomocí relativně malého počtu senzorů. 
    \item Identifikace permeability z kontinuálního měření tlaku je velmi citlivá a při větším množství tlakových senzorů umožňuje částečně zachytit i prostorovou variabilitu permeability.
    \item Citlivost na počáteční napětí, lze vhodně kombinovat s jinými metodami stanovení počátečního napětí a tím eliminovat možný negativní vliv velké prostorové variability.
    \item Kontinuální měření pevně osazenou aparaturou eliminuje nekonzistence při provádění časové posloupnosti měření. 
    \item Přestože senzory poskytují lokální informaci o~tlaku, dochází měřením přes komoru k~průměrování přes povrch komory a tím k~přirozenému potlačení šumu. 
\end{itemize}

\subsubsection{Slabé stránky návrhu}
\begin{itemize}
    \item Kontinuální měření pórovch tlaků lze smysluplně interpretovat pouze za přadpokladu plné saturace okolní horniny. Lokalita PVP Bukov je však ovlivněna dloudobým odčerpáváním  díky předešlé důlní činnosti. Kompaktní hornina bude vzhledem k malým hydraulickým vodivostem pravděpodobně stále saturována, avšak větší pukliny mohou již být jen částečně nasycené. V případě výskytu takových puklin může být interpretace měření výrazně komplikovanější. 
    \item Komplikovaná geologická stavba lokality patrně způsobuje značnou heterogenitu počátečního napětí. To může zhoršit podmínky pro citlivost inverzních metod.
    \item Větší počáteční investice: realizace měřících vrtů, jejich zapakrování a pečlivé osazení tlakových senzorů.
    \item Delší doba měření, omezení některých činností na lokalitě.
    \item Tlakové senzory nejsou při měření přístupné, a nelze tak opravit případné chyby, což zvyšuje nároky na kvalitu provedení.
\end{itemize}

\subsubsection{Příležitosti}
\begin{itemize}
\item Ověření modelů pro charakterizaci vlivu EDZ na indikátor bezpečnosti.
\item Možnost získání unikátních dat zahrnujících jak závislost hydraulické vodivosti na mechanických změnách, tak její prostorovou variabilitu. 
\item Při experimentu TSX byla pozorována vysoká korelace mezi hydraulickou vodivosti a mikroseismickým měřením vlnové rychlosti. Navrhovaný experiment umožní ověřit tuto korelaci v~českých podmínkách. Prokázání podobné korelace i na lokalitě úložiště by umožnilo identifikaci hydraulické vodivosti EDZ pomocí neinvazních a jednodušších seismických měření.
\item Ověření technologických schopností realizovat měření pro identifikaci EDZ v~dostatečné přesnosti.
\item Možné rozšíření či pokračování experimentu pro zahrnutí vlivu teploty na dlouhodobý vývoj EDZ.
\end{itemize}

\subsubsection{Rizika}
\begin{itemize}
    \item Malá citlivost měření v~případě příliš malého počátečního napětí, nebo příliš velké hydraulické vodivosti. V~tomto směru jsou podmínky na PVP Bukov středně příznivé.
    \item Velká prostorová variabilita počátečního napětí v~rámci PVP Bukov, bez prediktivního modelování zatím nejistý dopad na citlivost inverzních metod.
    \item Příliš komplikovaná geologická stavba a velká heterogenita vlastností horniny může být nad rozlišovací schopnost experimentu. Identifikovaná prostorová variabilita hydraulické vodivosti může pak být převážně nahodilá. Částečně lze eliminovat použitím pokročilých regularizačních technik při inverzní analýze.
    \item Větší množství tlakových senzorů a doprovodných měření zvyšuje riziko chyb. Pro snížení tohoto rizika je třeba dbát na kvalitu provedení (kabeláž, těsnění) a dokumentaci experimentu (značení kabelů, přesné informace o~polohách a časech měření, digitalizace získaných dat a pod.). Plán sběru a organizace získaných dat musí být součástí projektu experimentu.
    \item V~tuzemské komunitě jsou malé zkušenosti s propojením in-situ experimentů s komplexními numerickými simulacemi. Lze očekávat problémy v~komunikaci mezi více obory.
\end{itemize}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\bibliographystyle{plain}
%\bibliography{EDZ_plasticity.bib}
%\bibliography{references.bib}

% TRYING TO USE BIBLATEX
\printbibliography

\end{document}
