\section{Introduction}

We consider two-stage robust facility location problems under demand uncertainty where we are given a set of clients and a set of facilities in a common metric space. In the first stage, the decision-maker needs to select the (integral) units supply at each facility. The uncertain demand is then selected adversarially and needs to be satisfied by the existing supply with the minimum assignment cost in the second stage.  The goal is to determine the first-stage supply such that the sum of first-stage supply cost and the worst-case  assignment cost over all demand scenarios is minimized. Our problem is motivated by settings where the lead time to procure supply is large and obtaining additional units of supply in the second stage is not feasible. The uncertain second-stage demand must then be  satisfied by supply units from the first stage, a common constraint in many applications. This is a departure from existing work on two-stage robust facility location, network design and, more generally, robust covering problems that have been studied extensively in the literature~\cite{dhamdhere2005pay,FJMM07,anthony2010plant,gupta2014thresholded} where the second-stage decisions allow for ``adding more supply'' (more specifically adding more sets/facilities to satisfy the requirement).

In this paper, we consider  an  implicit  model  of  uncertainty  with exponentially-many demand scenarios specified by an upper bound $k$ on the number of second-stage demand clients.  Since the number of scenarios is exponentially-many (specified compactly), we can not efficiently solve even the LP relaxation for the problem. In contrast, if the set of second-stage scenarios are explicitly specified (for the explicit scenario model, see for instance \cite{dhamdhere2005pay,anthony2010plant}), we can write a polynomially-sized LP relaxation with assignment decisions for each scenario. The main challenge then is related to obtaining an integral solution, which for the case of set covering and several network design problems can be reduced to deterministic versions (see, for instance, Dhamdhere et al. \cite{dhamdhere2005pay}).

In contrast, in an implicit model of uncertainty (with possibly exponentially-many scenarios), one of the fundamental challenges is to even approximately solve the linear relaxation of the problem efficiently. The implicit model of uncertainty with an upper bound on number of uncertain second-stage clients or elements has been studied extensively in the literature.
 Feige et al. \cite{FJMM07} show under a reasonable complexity assumption that it is hard to  solve the LP relaxation of a two-stage set covering problem within a factor better than $\Omega( \log n / \log \log n)$ under this implicit model of uncertainty. They also give an $O(\log^2n)$-approximation for the $0-1$ two-stage robust set-covering problem. Gupta et al.~\cite{gupta2014thresholded} give an improved $O(\log n)$-approximation for the set-covering problem, thereby matching the deterministic approximation guarantee. El Housni and Goyal~\cite{housni2018optimality} show that a static policy (that is, linear in the set of second-stage elements, also referred to as an affine policy) gives an $O(\log n/\log\log n)$-approximation for the two-stage LP, thereby matching the hardness lower bound for the fractional problem. Gupta et al. \cite{gupta2016robust} and Khandekar et al. \cite{khandekar2008two} present approximations for several network design problems under the implicit model of uncertainty. Although there is a large body of work in this direction, as we mentioned earlier, our problem is different from set covering, since there is no possibility of adding recourse capacity;  prior results do not  imply an approximation for our model.

\noindent
{\bf Other related work.}
Many variants of  facility location problems have been studied extensively in the literature, including both deterministic versions as well as variants that address demand uncertainty. We refer the reader to \cite{shmoys2000approximation} for a  review of the many variants of deterministic facility location problems. Among the models that address demand uncertainty in the facility location problems, in addition to robust  \cite{atamturk2007,charikar2001algorithms}, there are also stochastic \cite{gupta2004boosted,immorlica2004costs,swamy2006approximation,ravi2004hedging} and distributionally robust \cite{linhares2019approximation,basciftci2019distributionally,delage2010distributionally} models that have been studied extensively in the literature. In a stochastic model, there is a distribution on the second-stage demand scenarios and the goal is to minimize the total expected cost.
A distributionally robust model can be thought of as a hybrid between stochastic and robust where the second-stage distribution is selected adversarially from a pre-specified set and the goal is to minimize the worst-case expected cost. We refer the reader to the survey \cite{snyder2006facility} for an extensive review of facility location problems under uncertainty.

\subsection{Our Contributions}

Let us begin with a formal problem definition. We are given a set of $n$ facilities $\cal F$ and $m$ clients $\cal C$ in a common metric $d$ where $d_{ij}$ denotes the distance between $i$ and $j$.
For each facility $i \in {\cal F}$, there is a cost $c_i$ per unit of supply at $i$.
The demand uncertainty is modeled by an implicit set of scenarios ${\cal C}^k$ that includes all subsets of clients $\cal C$ of size at most $k$.
The decision-maker needs to select an (integral) number of units of supply $x_i$ for each facility $i\in \cal F$ in the first-stage. An adversary observes the first stage decisions and selects a worst-case demand scenario $S \in {\cal C}^k$  that must be satisfied with the first-stage supply, where each client in the realized scenario needs one unit of supply. The goal is to minimize the sum of the first-stage supply cost and the worst-case  assignment cost over all second-stage demand scenarios. We refer to this problem as {\textit {soft-capacitated robust facility location}} \eqref{scrfl}. Typically in the literature, {\em soft-capacitated} refer to settings where violations of capacity upper bounds are allowed. The analogue here is that we can add any amount of supply in a facility without  upper bounds ($x_i \in \mathbb{Z_+})$ but we pay a per unit cost of supply. %

We consider a class of static assignment policies, where each of the $m$ clients has a static fractional assignment to facilities that is independent of the scenario, leading to a feasible second-stage solution for each demand scenario, while respecting supply capacities. Note this is a restriction, since the optimal second-stage assignment decisions are scenario-dependent in general. As a warm-up, we show that static assignment policies are optimal for the  uncapacitated case with unlimited supply at each open facility (i.e., there is a cost $c_i$ to open facility $i$ with unlimited supply). We refer to this problem as {\textit {uncapacitated robust facility location}} \eqref{urfl}. This is based on the intuition that each client can be assigned to the closest open facilities in an optimal solution in any scenario; this leads to  optimality of a static assignment policy for the LP relaxation (Theorem \ref{thm:static-opt}).

\begin{theorem} \label{thm:static-opt}
A static assignment policy is optimal for the linear relaxation of \eqref{urfl}.
\end{theorem}

The optimality of static assignment is not true in general when the supply at facilities is constrained (or equivalently, there is a cost per unit of supply). The main contribution in this paper is to show that static assignment policies give an $O(\log k/\log \log k)$-approximation for the LP relaxation of \eqref{scrfl} (Theorem~\ref{thm:static:log}). We show this by constructing such a solution, starting from an optimal first-stage supply. The optimal static assignment policies can be computed efficiently by solving a compact LP.

\begin{theorem}\label{thm:static:log}
A static assignment policy  gives $O(\log k / \log \log k)$-approximation for the linear relaxation of \eqref{scrfl}.
\end{theorem}

Furthermore, the fractional supply in the first stage can be rounded to an integral supply using ideas similar to rounding algorithms for the deterministic facility location \cite{shmoys1997approximation}. In particular, the static assignment solution for the uncapacitated case can be rounded to give a $4$-approximation algorithm for \eqref{urfl}.  The static assignment solution for the soft-capacitated  case can be rounded within a constant factor, which results in  an $O(\log k / \log \log k)$-approximation algorithm for \eqref{scrfl}. We would like to note that while the fractional assignment is static in our approximate LP solution, our integral assignment for any client in the second-stage  depends on the other demand clients in the scenario; thereby, making our static assignment policy adaptive in implementation.

\section{Warm-up: uncapacitated robust facility location  }

\subsection{Problem formulation}

In this section, we consider the {\textit{uncapacitated robust facility location problem}} ({\sf URFL}) where for each $i \in {\cal F}$, there is a cost $c_i$ to open facility $i$ with unlimited supply. The problem can be stated as the following integer program, where each binary variable $x_i, i \in {\cal F}$ indicates if  facility $i$ is opened and each $y_{ij}^S, i \in {\cal F}, j \in S, S \in {\cal C}^k$ indicates the assignment of client $j$  to facility $i$ in scenario $S$.

\begin{equation}\label{urfl} \tag{\sf URFL}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i +\max_{ S \in {\cal C}^k} \sum_{i \in {\cal F}} \sum_{ j \in {S}}  d_{ij} y_{ij}^S &\\
\text{s.t.    } & \sum_{i \in {\cal F}} y_{ij}^S \geq 1, \qquad &   \forall S \in {\cal C}^k,\forall j \in {S},     \\
& x_i \geq y_{ij}^S, \qquad  & \forall i  \in {\cal F},   \forall S \in {\cal C}^k,\forall j \in {S}, \\
& x_i \in \{0,1\}, \; y_{ij}^S \geq  0, \qquad & \forall i  \in {\cal F},   \forall S \in {\cal C}^k, \forall j \in {S}. \\
\end{aligned}
\end{equation}

Note that the second-stage problem is a transportation problem and since the demand of a client is integral (0 or 1), the  optimal solution $y^S_{ij}$ is integral as well.
The special case of \eqref{urfl} where the uncertainty set contains only a single scenario corresponds to the NP-hard classical and well-studied uncapacitated facility location problem, which is hard to approximate within a constant better than 1.463 unless NP has an  $O(n^{O(\log \log n) })$-time algorithm \cite{guha1999greedy}.
We let  ({\sf LP-URFL}) denote the linear relaxation of \eqref{urfl},  where we replace $x_i \in \{0,1\}$ by $x_i \geq 0$ for each $i \in {\cal F}$.
We would like to note that, in our model, the number of scenarios could be exponential in the dimension of the problem.
Hence, in general, even the linear relaxation of such a problem could be challenging. However, we show that  ({\sf LP-URFL})  can be solved in polynomial time using a  \textit{Static Assignment Policy} for the second-stage variables. Moreover, we can round the fractional solution losing only a constant factor, thereby  getting  a constant approximation for \eqref{urfl}. This section serves as a warm-up for introducing and motivating static assignment policies before addressing the class of capacitated robust facility location problems.

\subsection{Static assignment policy}

Consider an optimal solution of \eqref{urfl}. Since each open facility can have an unlimited amount of supply, each client in the realized scenario is assigned to the closest facility among the opened ones.  Therefore, a client $j$ is always assigned to the same open facility in all scenarios $S$ where $j \in S$.   The same observation holds as well for ({\sf LP-URFL}) where each client is assigned to the same fractionally opened facilities independent of the realized scenario. Thus, the assignment of a client is static. This can be captured by the following policy.

\vspace{2mm}
\noindent
{\bf Static assignment policy.} There exists $y_{ij} \geq 0$ for each $ i \in {\cal F}, j \in {\cal C}$ such that
\begin{equation}
 \label{eq:policy}
 \forall {S} \in {\cal C}^k,  \forall i \in {\cal F},\forall j \in S, \qquad \qquad  y_{ij}^S=    y_{ij}.
 \end{equation}

\noindent
{{\em Proof of Theorem \ref{thm:static-opt}}}. Let $(\mb x^*, \mb y^{*S}, S \in {\cal C}^k )$ be an optimal solution to  {\sf (LP-URFL)}. Since there are no capacities on facilities, each client $j$ is assigned to the closest fractionally opened facilities. In particular, for each $j \in {\cal C}$,
let $\pi_j$ be a permutation of ${\cal F}=\{1,\ldots,n\}$ such that  $ d_{\pi_j(1)j} \leq d_{\pi_j(2)j}  \leq \ldots \leq d_{\pi_j(n)j}$, and let
 $\ell = \min \{  p  \;   \vert  \;    x^*_{\pi_j(1)}+x^*_{\pi_j(2)}+ \ldots + x^*_{\pi_j(p)} \geq  1 \} $.
Denote $ \hat{x}_{\pi_j(\ell)} = 1- (x^*_{\pi_j(1)}+\ldots + x^*_{\pi_j({\ell-1)}}).$
The optimal solution can be written in the form \eqref{eq:policy} as follows: for $S \in {\cal C}^k$ and $j \in S$; ${y}_{ij}^S =   x^*_i$  for $i  \in \{\pi_j(1),\ldots,\pi_j(\ell-1) \}$, ${y}_{ ij}^S=  \hat{x}_i $  for $i  = \pi_j(\ell)$, and  ${y}_{ ij}^S=0 $, otherwise.
\qed

Let ({\sf Static-URFL}) denote the  problem after  restricting the second-stage variables $y_{ij}^S$ in {\sf (LP-URFL)} to a policy  \eqref{eq:policy}, which can then be reformulated as follows:
\begin{equation} \tag{{\sf Static-URFL}}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i  +\max_{ S \in {\cal C}^k} \sum_{i \in {\cal F}} \sum_{ j \in {\cal C}} {\bf 1} (j \in S) \cdot d_{ij} y_{ij} \\
 \text{s.t.    }\; & \sum_{i \in {\cal F}} y_{ij} \geq 1, \qquad & \forall j \in {\cal C}, \\
& x_i \geq y_{ij} \geq 0, \qquad  & \forall i  \in {\cal F},  \forall j \in {\cal C}. \\
\end{aligned}
\label{static-urfl}
\end{equation}

From Theorem \ref{thm:static-opt}, \eqref{static-urfl} is equivalent to {\sf (LP-URFL)}. The number of variables in \eqref{static-urfl} is reduced to a polynomial number since the $y_{ij}$ no longer depend  on the scenario $S$. The inner maximization problem is still taken over an exponential number of scenarios; however we can separate efficiently over these scenarios and write an efficient compact LP formulation for \eqref{static-urfl}:
\begin{equation} \label{eq:dual}
\begin{aligned}
 \max_{S \in {\cal C}^k } \; \{ \sum_{i \in {\cal F}} \sum_{j \in  {\cal C}} {\bf 1} (j \in S) \cdot d_{ij} y_{ij}    \}
    = \max_{\mb h \in [0,1]^{\vert {\cal C} \vert}} \; \{ \sum_{i \in {\cal F}} \sum_{j \in {\cal C}}  d_{ij} y_{ij} h_j  \;  \vert \;  \sum_{j \in {\cal C} } h_j \leq k \} \\
 =  \min_{\mu, \mb \omega \geq 0 } \; \{ k \mu + \sum_{j \in {\cal C}} \omega_j  \; \vert  \; \mu + \omega_j \geq  \sum_{i \in {\cal F}} d_{ij} y_{ij}, \; \forall j \in {\cal C} \},
\end{aligned}
\end{equation}
where the first equality holds because the optimal solution of the right maximization problem occurs at the extreme points of the $k$-ones polytope, which corresponds to the worst-case scenarios of ${\cal C}^k$ and the second equality follows from strong duality. Therefore, by dropping the min and introducing $\mu$ and all $  \omega_j$ as variables, we  reformulate \eqref{static-urfl} as the following  linear program:

\begin{equation} \label{eq:lp1}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i + k \mu + \sum_{j \in {\cal C}} \omega_j \\
 \text{s.t.    }\; & \mu + \omega_j \geq  \sum_{i \in {\cal F}} d_{ij} y_{ij}, \qquad & \forall j \in {\cal C}, \\
 &  \sum_{i \in {\cal F}} y_{ij}\geq 1, \qquad & \forall j \in {\cal C},     \\
& x_i \geq y_{ij},    \qquad  & \forall i  \in {\cal F},  \forall j \in {\cal C}, \\
& x_i \geq 0, \; y_{ij} \geq 0, \;  \omega_j \geq 0, \; \mu \geq 0,   \qquad  & \forall i  \in {\cal F},  \forall j \in {\cal C}.
\end{aligned}
\end{equation}

Finally, we  round  the solution of \eqref{static-urfl}   to an integral solution for \eqref{urfl} while losing  only a constant factor. This can be done using prior work on rounding techniques from the literature of deterministic facility location problems. In fact,  the LP rounding technique in Shmoys et al. \cite{shmoys1997approximation}, which gives a $4$-approximation algorithm to the deterministic uncapacitated problem also gives a $4$-approximation algorithm for \eqref{static-urfl}. The idea is to define a ball around each client of radius equal to the fractional assignment cost of the client (which is independent of any scenario for our static policy). Then, we open facilities in non-intersecting balls of ascending radius. The result is given in the following theorem and for completeness,  the details  of the rounding are in Appendix  \ref{appendix:A}.
\begin{theorem} \label{thm:rounding:uncapacitated}
\eqref{static-urfl}  can be rounded to give a 4-approximation to {\sf (URFL)}.
\end{theorem}
We would like to note that the focus in this section is not about finding the  best constant approximation  for {\sf (URFL)}, but we introduce it as a warm-up for motivating the static assignment policy before presenting our main result in the next section.

\section{Soft-capacitated robust facility location }

\subsection{Problem formulation}

In this section, we consider  the {\textit {soft-capacitated robust facility location}} ({\sf SCRFL}) which is similar to \eqref{urfl} except that each facility $i$ incurs a  linear supply cost,  where $c_i$ is the cost per unit of supply. We refer to $x_i$ as the supply (or capacity) in facility $i$. Each client in the realized scenario needs to be satisfied by one unit of supply. The problem is called soft-capacitated since there is no upper bound on  $x_i$, ($x_i \in \mathbb{Z_+})$. The problem is given by the following integer program:
\begin{equation}\label{scrfl} \tag{\sf SCRFL}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i +\max_{ S \in {\cal C}^k} \sum_{i \in {\cal F}} \sum_{ j \in {S}}  d_{ij} y_{ij}^S &\\
\text{s.t.    } & \sum_{i \in {\cal F}} y_{ij}^S \geq 1, \qquad & \forall S \in {\cal C}^k , \forall j \in {S},     \\
& x_i \geq \sum_{j \in S} y_{ij}^S,\qquad  & \forall i  \in {\cal F},   \forall S \in {\cal C}^k, \forall j \in {S}, \\
& x_i \in \mathbb{N}, \; y_{ij}^S \geq  0, \qquad & \forall i  \in {\cal F}, \forall S \in {\cal C}^k, \forall j \in {S} . \\
\end{aligned}
\end{equation}

We let  ({\sf LP-SCRFL}) denote the linear relaxation of \eqref{scrfl},  where we replace $x_i \in \mathbb{N}$ by $x_i \geq 0$, for each $i \in {\cal F}$. We would like to note that even the linear relaxation ({\sf LP-SCRFL}) is challenging to solve since it has exponentially-many variables (scenarios).  Unlike the uncapacitated case, the static assignment policy \eqref{eq:policy} is not  optimal for ({\sf LP-SCRFL}) and the optimal assignment for each client depends, in general, on the realized scenario. In particular, the same client could be assigned to  different facilities in  different scenarios. In contrast, we show  the surprising result that a static assignment policy gives $O( \log k/ \log \log k)$-approximation  to ({\sf LP-SCRFL}). Moreover,  we can round the solution of the static assignment policy to an integral solution for \eqref{scrfl} and  only lose an additional constant factor.
We  let ({\sf Static-SCRFL}) denote the problem when we  restrict the second-stage variables $y_{ij}^S$ in ({\sf LP-SCRFL}) to static assignment policies  \eqref{eq:policy}. The problem  can then be reformulated as follows:

\begin{equation}\label{static-scrfl} \tag{\sf Static-SCRFL}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i +\max_{ S \in {\cal C}^k} \sum_{i \in {\cal F}} \sum_{ j \in {\cal C}} {\bf 1} (j \in S) \cdot  d_{ij} y_{ij} &\\
\text{s.t.    } & \sum_{i \in {\cal F}} y_{ij} \geq 1, \qquad & \forall j \in {\cal C},  \\
& x_i \geq  \max_{ S \in {\cal C}^k}     \; \;  \sum_{ j \in  {\cal C}} {\bf 1} (j \in S)    \cdot y_{ij}, \qquad & \forall i \in {\cal F}, \\
& x_i \geq 0, \; y_{ij} \geq  0, \qquad & \forall i  \in {\cal F},  \forall j \in {\cal C}. \\
\end{aligned}
\end{equation}

\subsection{An $O ( \frac{ \log k} { \log \log k} )$-approximation algorithm}

Our main contribution in this section is to show that a static assignment policy \eqref{eq:policy}  gives $O(\log k / \log \log k)$-approximation  for {\sf (LP-SCRFL)} (Theorem~\ref{thm:static:log}).
To prove this theorem, we consider an optimal solution of ({\sf LP-SCRFL}) and massage it to construct a solution of the form  \eqref{eq:policy} while losing  $O( \log k / \log \log k)$ factor. We first present our construction and several structural lemmas and then give the proof of Theorem \ref{thm:static:log}.

\vspace{2mm}
\noindent
{\bf Our construction.}
Let $\mb {x^*}: (x^*_i)_{ i \in {\cal F}}$ be an optimal first-stage solution of ({\sf LP-SCRFL}), let $\mathrm{OPT}_1$ be the corresponding optimal first-stage cost and let $\mathrm{OPT}_2$ be the corresponding optimal second-stage cost. We will classify the clients ${\cal C}$ into three subsets $C_1,C_2,C_3$ using Procedure \ref{alg-soft} (below) and then specify a static assignment policy for each subset. We use the following notation in the procedure.
Let  $\alpha > 1$ and $r= 5 \cdot \mathrm{OPT}_2/ k$. For $\ell \geq 1$ and $j \in {\cal C}$, we let $B_j^{\ell}$ denote the ball centered at client $j$ of radius $\ell r$. We initialize the sets $F \leftarrow {\cal F}$ and $C \leftarrow {\cal C}$ and update them at iteration, as explained in the procedure, until $C$ becomes empty. We let $Cl(B)$ denote the set of clients in $C$ that are inside the ball $B$ and let $Sp(B)$ denote  the total optimal supply of facilities $F$ that are inside the ball $B$, i.e.,
$$ Sp(B) = \sum_{i \in  { F}} {\bf 1} (i \in B) \cdot x_i^* \qquad  \text{and} \qquad Cl(B) =  \{ j \in  { C} \; \vert \; j \in   B \}. $$
Note that both $ Sp(B)$ and $Cl(B) $ depend on the current sets of facilities $F$ and clients $C$, which we update at each iteration of the while loop in the procedure. But for the ease of notation, we do not refer to  them  with the indices $F$ and $C$.

In the procedure, while the set $C$ is not empty, we pick a client $j \in C$ and grow three balls around it: $B_j^{2 \ell - 1}$ ({\em internal} ball), $B_j^{2 \ell }$ ({\em medium} ball) and $B_j^{2 \ell + 1}$ ({\em external} ball) starting with $\ell=1$. For each $\ell$, we check if the number of clients in the internal ball $B_j^{2 \ell - 1}$ is greater than $k$ (line 4);  if this is the case, we remove them from $C$, put them in $C_1$ and restart in line 2. If not, we check if the supply in the medium ball $B_j^{2 \ell}$ is  sufficient to satisfy half of the clients in the internal ball $B_j^{2 \ell - 1}$ (line 7); if that is not the case, we remove those clients from $C$, put them in $C_2$, and restart in line 2. Otherwise, we finally check if the supply in the the medium ball $B_j^{2 \ell}$ is sufficient to satisfy a fraction $1/2 \alpha$ of the clients in the external ball  $B_j^{2 \ell+1}$ (line 10); if that is the case we remove all the clients in $B_j^{2 \ell+1}$ and put them in $C_3$, we also  remove all the facilities in $B_j^{2 \ell}$ and restart in line 2. If none of these three conditions holds, we increase $\ell$ to $\ell+1$. First, we show that after at most $\log_{\alpha} k$ increments (i.e., $\ell \leq \log_{\alpha}k $), one of three conditions must hold and therefore we will remove some clients from $C$ and restart in line 2. Which implies that after a finite number of iterations, the set $C$ becomes empty. In particular, We have the following lemma.

\floatname{algorithm}{Procedure}

\begin{algorithm*}[t]
\caption{}\label{alg-soft}
\begin{algorithmic}[1]
\State  Initialize  $C \leftarrow {\cal C}, C_1 \leftarrow \emptyset, C_2 \leftarrow \emptyset, C_3 \leftarrow \emptyset, F \leftarrow {\cal F} $.
\While{$ { C} \neq \emptyset$ }
\State Pick a client $j \in C$. Initialize $\ell=1$
\If {  $  \vert   Cl(B_j^{2 \ell -1}) \vert \geq k$    }
\State $C_1  \leftarrow C_1 \cup Cl(B_j^{2 \ell -1})$, \; $ { C} \leftarrow  { C} \setminus Cl(B_j^{2 \ell -1})$
\State Stop, return to line 2
\EndIf
\If {  $  Sp(B_j^{2 \ell })         <  \frac{1}{2}  \cdot \vert   Cl(B_j^{2 \ell - 1})  \vert $    }
\State $C_2 \leftarrow C_2 \cup Cl(B_j^{2 \ell - 1})$, \;  $ { C} \leftarrow   {C} \setminus Cl(B_j^{2 \ell - 1})$
\State Stop, return to line 2
\EndIf
\If {  $  Sp(B_j^{2 \ell })         \geq   \frac{1}{2 \alpha}   \cdot     \vert   Cl(B_j^{2 \ell +1})  \vert $    }
\State $C_3 \leftarrow C_3 \cup Cl(B_j^{2 \ell +1})$, \; $ { C} \leftarrow  { C} \setminus Cl(B_j^{2 \ell +1})$,
\State$ { F} \leftarrow { F} \setminus \{ i \in { F} \; \vert \; i \in  B_j^{2 \ell } \}$
\State Stop, return to line 2
\Else \;  {$ \ell \leftarrow \ell +1$, return to line 4.}
\EndIf
\EndWhile
\end{algorithmic}
\end{algorithm*}

\begin{lemma} \label{lem:steps}
In Procedure   \ref{alg-soft},  after a finite number of iterations, the set $C$  becomes empty and $C_1 \cup C_2 \cup C_3$ is equal to $ {\cal C}$. Moreover,   $\ell$ is always less than $  \log_{\alpha} k $.
\end{lemma}

\begin{proof}
Fix a client $j$ and let $\ell \geq 1$. If none of the three conditions (``if'' statements) holds then
$$     \alpha \cdot \vert   Cl(B_j^{2 \ell -1}) \vert   \leq 2 \alpha \cdot Sp(B_j^{2 \ell })         <          \vert   Cl(B_j^{2 \ell +1})  \vert . $$
Therefore, the number of clients grows geometrically when we increase the radius of the balls and by induction, we have that
$$ \alpha^{\ell}  \leq  \alpha^{\ell} \cdot \vert   Cl(B_j^{1}) \vert   <         \vert   Cl(B_j^{2 \ell +1}) \vert, $$
where $\vert   Cl(B_j^{1}) \vert \geq 1$, since  $Cl(B_j^{1})$  contains at least the client $j$.
Hence, after at most $\log_{    \alpha} k$ increments, we will reach $k$ clients,  and must stop by the first condition, and return to line 2. Hence, we always have $\ell \leq \log_{\alpha} k $. Finally since, we remove at least one client at each iteration of the while loop,  the set $C$ becomes empty after at most $\vert {\cal C} \vert$ iterations, and finally $C_1 \cup C_2 \cup C_3= {\cal C}$.
\qed
\end{proof}

Now we are ready to present our static assignment policy for {\sf(LP-SCRFL)}. The following three lemmas show our constructed static assignment for each client in the three subsets $C_1,C_2,C_3$. Moreover, we specify the supply used to satisfy each subset of these clients and present the analysis for the  assignment cost.

For each client in the set $C_1$, we know that it belongs to a ball with at least $k$ clients. By the feasibility of the optimal solution, this implies that there exists $k$ units of supply in $\mb x^*$ ``close'' to this ball. Hence, we  satisfy this client by using the same fraction $1/k$ of  $\mb x^*$ (static assignment)  while paying a small assignment cost (roughly a constant times the radius of the ball). Since, there are at most $k$ clients in each scenario, and each one is using at most $\mb x^*/k$, we  need to dedicate only one $\mb x^*$ for all clients in $C_1$.  Formally, we have the following lemma.

\begin{lemma} There exists a static assignment policy for $C_1$ such that each client in $C_1$ is using at most the supply $\mb {x^*}/k$  and has an assignment cost less than  $O(\log_{\alpha}k/k ) \cdot \mathrm{OPT}_2 $, i.e., there exists $(\tilde{y}_{ij})_{i \in {\cal F}, j \in { C_1}}$ such that  for each $j \in C_1$ :
$$         \sum_{i \in {\cal F}} \tilde{y}_{ij} \geq 1,  \qquad   \frac{x_i^*}{k}   \geq  \tilde{y}_{ij} \geq 0 , \; \forall i \in {\cal F},
\; \; \; \;\text{and} \; \; \; \; \sum_{i \in {\cal F}} d_{ij} \tilde{y}_{ij} = O(\log_{\alpha}k ) \cdot \frac{\mathrm{OPT}_2}{k}.$$
\label{lem:C1}
\end{lemma}

\begin{proof}
Let $j$ be a client of $C_1$. It is sufficient to show that the following minimization problem is feasible and its optimal cost is  $O(\log_{\alpha}k/k ) \cdot \mathrm{OPT}_2 $. Consider

\begin{equation} \label{eq:ss}
\min \; \left\{ \sum_{i \in {\cal F}} d_{ij} y_{ij} \; \bigg\vert \;        \sum_{i \in {\cal F}} y_{ij}  \geq  1, \; \; \frac{x_i^*}{k}       \geq  y_{ij} \geq 0 ,\; \; \forall i \in {\cal F} \right\} .
\end{equation}
Problem \eqref{eq:ss} must be feasible since the total supply in $\mb x^*$ is greater than the total demand in any scenario, i.e., $ \sum_{i \in {\cal F}} x_i^* \geq k$.
Recall that a client $j$ in $C_1$  belongs to one of the sets $Cl(B_{t}^{2 \ell - 1})$ for some  $t \in{\cal C}$  and $\ell \leq \log_{\alpha} k $ (Lemma \ref{lem:steps}) such that $\vert Cl(B_{t}^{2 \ell - 1}) \vert \geq k$. Consider a scenario $S$ formed by $k$ clients from $Cl(B_{t}^{2 \ell - 1})$. Let denote $\mb y^S$ the assignment of scenario $S$ in the optimal solution. Consider the following candidate solution for \eqref{eq:ss}:

$$ y_{ij}= \frac{1}{k} \cdot \sum_{ p \in S} y^S_{ip}, \; \forall i \in{\cal F}.$$
We have, by the feasibility of the optimal solution, for each $i \in {\cal F}$,
$ 0 \leq y_{ij} \leq  \frac{1}{k} x_i^* $
and
$$\sum_{i \in {\cal F}} y_{ij}  =   \frac{1}{k} \cdot       \sum_{i \in {\cal F}} \sum_{ p \in S} y^S_{ip}  \geq   \frac{1}{k} \sum_{ p \in S} 1 = 1.   $$
Therefore, our solution is feasible for \eqref{eq:ss}. Moreover,  we have
\begin{align*}
 \sum_{i \in {\cal F}} d_{ij}  y_{ij}  &= \frac{1}{k} \cdot \sum_{i \in {\cal F}} \sum_{ p \in S} d_{ij} y^S_{ip} \\
 & \leq   \frac{1}{k} \cdot  \sum_{i \in {\cal F}}  \sum_{ p \in S}  d_{ip} y^S_{ip}+ \frac{1}{k} \cdot \sum_{ p \in S}  d_{pj}         \\
 & \leq \frac{1}{k} \cdot \mathrm{OPT}_2   + 2(2\ell -1) r \\
 & \leq \frac{\mathrm{OPT}_2}{k}  + 2(2 \log_{\alpha}k      -1) \cdot 5 \cdot \frac{\mathrm{OPT}_2}{k}  = O(\log_{\alpha}k ) \cdot \frac{\mathrm{OPT}_2}{k},
\end{align*}
where the first inequality follows from the triangle inequality and the fact that
$\sum_{i \in {\cal F}}  y_{ip}^S=1$ for all $p \in S$ in the optimal solution. For the second inequality, we use the definition of $\mathrm{OPT}_2$ to bound the first term,  the second term $d_{pj}$ is bounded by the diameter of the ball $B_j^{2 \ell -1}$ which contains client $j$ and  all clients $p \in S$.
\qed
\end{proof}

Now, consider the set $C_2$. By construction, these are clients such that there is not  enough  supply  within a distance $r= 5\mathrm{OPT}_2/k$ to satisfy half of them. Therefore, intuitively they need to pay ``large'' distances  in the optimal assignment cost if they show up all together in the same scenario. In the following lemma, we show that we can  have no more than $k$ of these clients. As we would show later, this would imply that we can dedicate a supply $\mb x^*$ to $C_2$ and make a static assignment of all the clients $C_2$ to this $\mb x^*$.
\begin{lemma} The set $C_2$ has at most  $k$ clients.
\label{lem:C2}
\end{lemma}

\begin{proof}
Suppose, for the sake of contradiction, that $\vert C_2 \vert >  k$. Let $ G_1, G_2, \ldots, G_T$ be the disjoint  subsets of clients added at each iteration in the construction of $C_2$ in Procedure \ref{alg-soft}. In particular,  $C_2= G_1 \cup G_2 \cup \ldots \cup G_T$ for some $T$ where:
\begin{itemize}
    \item [(i)] for $t=1,2,\ldots,T, $ $G_t=Cl(B_{j_t}^{2 \ell_t-1})$ for some  client $j_t$ and  $  1  \leq \ell_t  \leq \log_{\alpha} k $.
    \item  [(ii)] the supply $Sp(B_{j_t}^{2 \ell_t})$ is less than half of the clients in $G_t$.
    \item [(iii)] each set $G_t$ has strictly less than $k$ clients, since the procedure has to fail the first ``if'' statement before adding $G_t$ into $C_2$.
\end{itemize}
Recall that $$ Sp(B_{j_t}^{2 \ell_t }) = \sum_{i \in  { F}} {\bf 1} (i \in B_{j_t}^{2 \ell_t }) \cdot x_i^*, $$
where $F$ is the current set of facilities in the procedure (and is not all $\cal F$ since some facilities have been removed in line 12 of the procedure). However, we would like to emphasize that when a  facility has been removed (in line 12 of the procedure),  all clients within distance  $r$ from this facility were removed as well (line 11). This is true, since  when we remove the facilities in a medium ball, say $B_j^{2 \ell }$, (line 12), we remove all  clients in the corresponding external ball $B_j^{2 \ell+1 }$ (line 11). Hence, the remaining clients in $C$ are at least $(2\ell+1)r - ( 2 \ell)r =r$ away from the removed facilities. In particular, for the clients $G_t$, the supply that has been removed before they were added to $C_2$ is at least $r$ away from them. Therefore, all the facility in $\cal F$ that are within a distance $r$ from a client in $G_t$ belong to the set $F$ that verifies $\sum_{i \in  { F}} {\bf 1} (i \in B_{j_t}^{2 \ell_t }) \cdot x_i^* \leq 2 \cdot \vert G_t \vert $. This implies that the  supply of all facilities within a distance $r$ from $G_t$ in the optimal solution, is less than half of the clients $G_t$. Hence, if all of the clients $G_t$ show up in a scenario,  the optimal second-stage solution needs to pay an assignment cost of at least $r \cdot \vert G_t \vert /2$.

Order the sets $G_t$ according to their cardinalities: wlog assume that $\vert G_1 \vert \geq \vert G_2 \vert \geq \ldots \geq \vert G_T \vert$.
We construct a scenario $\hat{S}$ by taking clients from the sets $G_1$, $G_2, \ldots$ until we hit $k$. This is possible since by assumption $\vert C_2 \vert >  k$. Assume that
$$ \vert G_1 \vert + \vert G_2 \vert + \ldots \vert G_{p-1}  \vert    +  \vert \bar{G}_{p}  \vert = k ,$$
for some $p$, where $ 2 \leq p \leq T$. Note that $\bar{G}_{p}$ is a subset of ${G}_{p}$, since we can reach $k$ before taking all the clients of the last set ${G}_{p}$.
For each $t=1, \ldots, p-1$, the optimal second-stage decision needs to pay at least $r \cdot \vert G_t \vert /2$. Therefore,

  $$ \mathrm{OPT}_2 \geq  \frac{1}{2} \cdot  r \cdot ( \vert G_1 \vert + \vert G_2 \vert + \ldots \vert G_{p-1}  \vert   ) .$$
We did not include  $  G_{p} $, since not all these clients  are necessary in the scenario $\hat{S}$, but only $\vert \bar{G}_{p}  \vert$ of them. Since $\bar{G}_{p}$ has the smallest cardinality
$$  \vert G_1 \vert + \vert G_2 \vert + \ldots \vert G_{p-1}  \vert \geq     \frac{1}{2} ( \vert G_1 \vert + \vert G_2 \vert + \ldots \vert G_{p-1}  \vert    +  \vert \bar{G}_{p}  \vert )  = \frac{k}{2}. $$
Therefore,
$$  \mathrm{OPT}_2 \geq      r \cdot \frac{k}{4}   =   5 \cdot \frac{\mathrm{OPT}_2}{4}  ,   $$
  which is a contradiction.
 Therefore, $\vert C_2 \vert \leq k$.
 \qed
\end{proof}

Finally, for clients $C_3$, we show that there exists $ \vert C_3 \vert /2\alpha$ units of supply ``close'' to them. In particular, we can  multiply these units by $2 \alpha$, dedicate them to $C_3$ and make a static assignment for $C_3$. We have the following lemma.

\begin{lemma}
There exists a supply $\mb {\hat{x}}$  that has a cost at most  $2 \alpha \cdot \mathrm{OPT}_1$ and there exists a static assignment policy such that all clients in $C_3$  are assigned to  supply $\mb {\hat{x}}$ and each client in $C_3$ has an assignment cost that is  $O(\log_{\alpha}k /k ) \cdot \mathrm{OPT}_2$, i.e., there exists $(\hat{y}_{ij})_{i \in {\cal F}, j \in { C_3}}$ and $(\hat{x}_i)_{i \in{\cal F}}$
 such that  for all $j \in C_3$ :
$$         \sum_{i \in {\cal F}} \hat{y}_{ij} \geq 1,   \qquad  \hat{x}_i  \geq \sum_{j \in C_3} \hat{y}_{ij},\; \forall i \in {\cal F},   \qquad    \hat{x}_i \geq 0, \; \hat{y}_{ij} \geq 0 , \; \forall i \in {\cal F},$$
$$\sum_{i \in {\cal F}} c_i   \hat{x}_i \leq 2 \alpha  \cdot \mathrm{OPT}_1 \qquad
\text{and} \qquad  \sum_{i \in {\cal F}} d_{ij} \hat{y}_{ij} = O(\log_{\alpha}k ) \cdot \frac{\mathrm{OPT}_2}{k}.
$$
 \label{lem:C3}
\end{lemma}

\begin{proof}
Let $ G_1, G_2, \ldots, G_T$ be the disjoint  subsets of clients added at each iteration to  construct  $C_3$ in Procedure \ref{alg-soft}. In particular,  $C_3= G_1 \cup G_2 \cup \ldots \cup G_T$ for some $T$ such that:  for all $t=1,2,\ldots,T, $ $G_t=Cl(B_{j_t}^{2 \ell_t+1})$ for some  client $j_t$ and some integer $\ell_t$ with $  1  \leq \ell_t  \leq \log_{\alpha} k $.  Moreover, the supply $Sp(B_{j_t}^{2 \ell_t})$ is greater  than a $1/2\alpha$ fraction of the clients in $G_t$. Hence, for each ball $B_{j_t}^{2 \ell_t}$, we multiply the supply by $2 \alpha $, move it to the cheapest facility in this ball and make a static assignment of all clients  $G_t$ to this cheapest facility. Since the supply in $B_{j_t}^{2 \ell_t}$ is removed along with clients $G_t$, it will not be used by the other clients in $C_3$.

Formally, let $i_t$ be the cheapest facility in the ball $B_{j_t}^{2 \ell_t}$. We define, for each facility $i$ in $B_{j_t}^{2 \ell_t}$, $\hat{x}_{i} = 2 \alpha \sum_{i' \in  { F}} {\bf 1} (i' \in B_{j_t}^{2 \ell_t}) \cdot x_{i'}^* $ if $i =i_t$ and $\hat{x}_{i}=0$ otherwise. For each client $j \in C_3$, we let $\hat{y}_{ij}=1$ for $i=i_t$ and $j \in G_t$, and let $\hat{y}_{ij}=0$, otherwise. Therefore, the first desired constraints in the lemma are verified. Let us check the last one. The distance between a client and its assigned facility in our solution is at most $r$ plus the diameter of the ball $B_{j_t}^{2 \ell_t}$, i.e., $$r+ 4 {\ell}_t r \leq (4 \log_{\alpha} k +1) \cdot 5 \cdot \mathrm{OPT}_2/k = O( \log_{\alpha}k/k) \cdot \mathrm{OPT}_2.$$
\qed
\end{proof}

\noindent
{\em Proof of Theorem \ref{thm:static:log}.}
Let $(\tilde{y}_{ij})_{i \in {\cal F}, j \in { C_1}}$ be the solution given in Lemma \ref{lem:C1} for satisfying the clients in $C_1$. We dedicate a supply $\mb x^*$ to clients $C_1$. Let $(\hat{y}_{ij})_{i \in {\cal F}, j \in { C_3}}$ and $(\hat{x}_i)_{i \in{\cal F}}$ be the solution given in Lemma \ref{lem:C3} for satisfying the clients in $C_3$. Finally, we know from Lemma \ref{lem:C2} that $C_2$ has at most $k$ clients, and therefore $C_2$ is a scenario. So we dedicate a supply $\mb x^*$ to $C_2$ and let the optimal assignment $y_{ij}^{C_2}$ be  our static assignment solution for $C_2$. In particular, we give the following solution to {\sf (LP-SCRFL)}, where the first stage solution is $  2 \mb x^* +   \mb {\hat{x}} $ and the static assignment policy is for all $ i \in {\cal F}$: $y_{ij} = \tilde{y}_{ij} $  for $ j \in C_1 $, $ y_{ij}= {y}^{C_2}_{ij}$ for $j \in C_2 $, and  $ y_{ij}=\hat{y}_{ij} $ for $ j \in C_3$. It is clear that $ \sum_{i \in {\cal F}} y_{ij} \geq 1 $, for each $j$ in $C_1 \cup C_2 \cup C_3$. Moreover, for any scenario $S \in {\cal C }^k$ and $i \in {\cal F}$,
\begin{align*}
    \sum_{j \in S} y_{ij} &=  \sum_{j \in S \cap C_1} \tilde{y}_{ij} + \sum_{j \in S \cap C_2} y^{C_2}_{ij} + \sum_{j \in S \cap C_3} \hat{y}_{ij} \\
    &\leq  \sum_{j \in S\cap C_1} \frac{x_i^*}{k} + x_i^* + \hat{x}_i \leq  2x_i^* + \hat{x}_i.
\end{align*}
Therefore, our solution is feasible for {\sf (LP-SCRFL)}. Let us evaluate its cost. The cost of the first stage is at most  $2 OPT_1 + 2 \alpha \mathrm{OPT}_1= O(\alpha) \cdot\mathrm{OPT}_1$.
For the second-stage cost, consider any scenario $S \in {\cal C}^k$,
We have
\begin{align*}
    \sum_{i \in {\cal F}} \sum_{ j \in S} d_{ij} y_{ij} &=  \sum_{j \in S \cap C_1}  \sum_{i \in {\cal F}} d_{ij}  \tilde{y}_{ij} +  \sum_{j \in S \cap C_2}   \sum_{i \in {\cal F}} d_{ij}y^{C_2}_{ij} + \sum_{j \in S \cap C_3} \sum_{i \in {\cal F}} d_{ij}\hat{y}_{ij} \\
    &\leq  \sum_{j \in S\cap C_1} O(\log_{\alpha}k ) \cdot \frac{\mathrm{OPT}_2}{k} + \mathrm{OPT}_2 + \sum_{j \in S \cap C_3} O(\log_{\alpha}k ) \cdot \frac{\mathrm{OPT}_2}{k} \\
    &\leq  O(\log_{\alpha}k ) \cdot \mathrm{OPT}_2 + \mathrm{OPT}_2 +  O(\log_{\alpha}k ) \cdot \mathrm{OPT}_2  \\
    &= O(\log_{\alpha}k ) \cdot \mathrm{OPT}_2.
\end{align*}
By balancing the terms $\alpha $ and $\log_{\alpha} k$, we choose $ \alpha =\log k / \log \log k$ which gives $O(\log k / \log \log k)$-approximation to {\sf (LP-SCRFL)}.
\qed

Similar to the uncapacitated problem, we can solve \eqref{static-scrfl} efficiently using a compact linear program. In fact, we dualize the inner maximization problem in the objective function of \eqref{static-scrfl} in the same way as \eqref{eq:dual}. In addition to that, we reformulate the second constraint in \eqref{static-scrfl} using the same dualization technique as follows:  for each $i \in {\cal F}$,

\begin{align*}
 \max_{S \in {\cal C}^k } \; \{ \sum_{ {j \in  S}}   y_{ij}  \}
 &   = \max_{\mb h \in [0,1]^{\vert C \vert}} \; \{ \sum_{ {j \in {\cal C}}}   y_{ij} h_j  \; \; \vert \; \; \sum_{j \in {\cal C} } h_j \leq k \} \\
 &=  \min_{\eta_i, \lambda_{ij} \geq 0 } \; \{ k \eta_i + \sum_{j \in {\cal C}} \lambda_{ij}  \; \vert  \; \eta_i + \lambda_{ij} \geq    y_{ij}, \; \forall j \in {\cal C} \}.
\end{align*}
The linear program is given by
\begin{equation}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i + k \mu + \sum_{j \in {\cal C}} \omega_j \\
 \text{s.t.    }\; & \mu + \omega_j \geq  \sum_{i \in {\cal F}} d_{ij} y_{ij}, \qquad & \forall j \in {\cal C}, \\
 &  \sum_{i \in {\cal F}} y_{ij}\geq 1, \qquad & \forall j \in {\cal C} ,    \\
  & x_i \geq k \eta_i + \sum_{j \in {\cal C}} \lambda_{ij} ,& \forall i \in {\cal F} ,    \\
    &\eta_i +  \lambda_{ij} \geq y_{ij}, & \forall i \in {\cal F} , \forall j \in {\cal C},   \\
& x_i, y_{ij}, \lambda_{ij}, \eta_i, \omega_j, \mu \geq 0,  \qquad  & \forall i  \in {\cal F},  \forall j \in {\cal C}, \\
\end{aligned}
\end{equation}
which can be reduced, after removing the variables $y_{ij}$, to
\begin{equation} \label{eq:lp2}
\begin{aligned}
\min \; &  \sum_{i \in {\cal F}} c_i x_i + k \mu + \sum_{j \in {\cal C}} \omega_j \\
 \text{s.t.    }\; & \mu + \omega_j \geq  \sum_{i \in {\cal F}} d_{ij} (\eta_i +  \lambda_{ij}) , \qquad & \forall j \in {\cal C}, \\
 &  \sum_{i \in {\cal F}} \eta_i +  \lambda_{ij}\geq 1, \qquad & \forall j \in {\cal C} ,    \\
  & x_i \geq k \eta_i + \sum_{j \in {\cal C}} \lambda_{ij}, & \forall i \in {\cal F}  ,   \\
& x_i, \lambda_{ij}, \eta_i, \omega_j, \mu \geq 0,  \qquad  & \forall i  \in {\cal F},  \forall j \in {\cal C}. \\
\end{aligned}
\end{equation}

Finally, we  round the optimal solution of \eqref{static-scrfl}  to an integral solution using the filtering and rounding techniques from Shmoys et al. \cite{shmoys1997approximation} while losing only a  factor of 12. Again, this rounding technique was designed for the deterministic facility location problem, but the same argument works as well for \eqref{static-scrfl}. Finally, since  \eqref{static-scrfl} gives $O(\log k / \log \log k)$-approximation to ({\sf LP-SCRFL}) and we only loose a constant factor in the rounding, this results in $O(\log k / \log \log k)$- approximation algorithm for ({\sf SCRFL}). We state the result in the following theorem and, for completeness, we present  the details  of the rounding in Appendix \ref{appendix:B}.
\begin{theorem} \label{thm:rounding:capacitated}
\eqref{static-scrfl} can be rounded to  give  $O(\frac{\log k}{\log \log k})$-approximation algorithm to \eqref{scrfl}.
\end{theorem}
Note that after rounding the supply in our solution of \eqref{static-scrfl},  the integral second-stage assignment for each realized scenario is a transportation problem and therefore its optimal solution is integral. We would like to emphasize that while the fractional assignment in our solution is static,  our integral assignment is not  necessarily static. In fact, a policy with an integral static assignment could even be bad our model.

\section{Conclusion.}
In this paper, we give a $O(\log k / \log \log k)$-approximation for soft-capacitated robust facility location problems with an implicit model of demand uncertainty. It is an interesting open question to study whether there exists  a constant approximation algorithm for the problem, even in special cases such as the Euclidean metric. Our solution approach relies on static fractional assignment policies, which we show are optimal for the  uncapacitated problem and give a strong theoretical guarantee for soft-capacitated case.
Static assignment policies, while reasonable for the case of soft-capacities can be shown to be arbitrarily bad for the case of hard-capacities, where in addition to cost per unit, there is also an upper bound on supply at each facility. It is another interesting open direction to study any non-trivial approximation in this setting.

{\small
{\newpage
\bibliographystyle{abbrv}


\section{Introduction}

Faults in the Earth's crust are complex systems along which earthquakes nucleate and propagate (e.g., \cite{Wibberley_complex}). Faults hold structures and heterogeneities at all scales (\cite{fault_imaging, fault_imaging2, rock_surf_bandwidth}). While they are often simplified to their simplest two-dimensional description (i.e., the fault plane), increasing complexity is now added to faults models (e.g., \cite{Rice_complex}). It is indeed considered that, to fully understand seismicity in various areas, it is paramount to account for some disorder in the faults frictional properties such as secondary faulting, off-fault damage or roughness of the fault plane (\cite{hetero_1, hetero_2, hetero_3, hetero_4, interlocked_rough_faults, fault_imaging_3, roughness1, roughness2}). For instance, the volume of damaged rocks, during the activation of a fault, depends on the initial contact roughness (\cite{wear_volume}) and, thus, a fault with a stronger roughness presents a different energy budget than a flat fault, as more energy is converted into surface area energy. In particular, roughness encourages the triggering of local events, but is believed to prevent the propagation of large-slip earthquakes (\cite{roughness2}). Additionally, large scale roughness tends to inhibit the propagation of any rupture faster than the shear wave velocity of surrounding rocks (\cite{bouchon_faulting_2010}).\\
Another degree of complexity is more rarely considered when modelling geological contacts and fault slip: the possible anisotropy in their frictional properties.
Morphological anisotropy is a known feature of faults, notably impacting the seismic waves velocity in their vicinity (\cite{seismic_anisotropy_3, seismic_anisotropy_1, seismic_anisotropy_2}) or the mobility of natural and injected fluids (\cite{perm_anisotropy}) in the subsurface. Frictional anisotropy, interestingly, is also regularly studied in other fields than seismology, for instance the tribology of rubber tires (\cite{rubber_1, rubber_2}), the strength of advanced adhesives (\cite{aniso_adhesive}), or the mitigation of water condensation (\cite{Condensation}). It is also considered to play a major role in nature (\cite{bio_friction_aniso}), for instance in the motion of numerous animals (\cite{aniso_adhesive, snake_friction_aniso, butterfly}) and the hydration of some plants (\cite{grass, lotus_bamboo}). In most cases, frictional anisotropy derives from the existence of preferential topographical orientations on, at least, one of the contact surfaces (\cite{friction_aniso_general, friction_anisotropy_2}). The length scale for such structural directivity can be as small as micrometer (\cite{friction_aniso_micro}) to nanometer (\cite{friction_aniso_atom, friction_aniso_atom_2}).
\begin{figure}
  \includegraphics[width=1\linewidth]{corona_metric.png}
  \caption{Topography (i.e. roughness) of the Corona Heights fault at the meter scale (\cite{fault_database, doi_fault_database}). This surface, called $S_{m}$, has a radius of $1$\,m and is defined on a $5$\,mm grid with a $1.25$\,mm elevation resolution. A parametric angle $\theta$ is defined from the main groove orientation.}
  \label{corona_metric}
\end{figure}
\begin{figure}
  \includegraphics[width=1\linewidth]{corona_millimetric.png}
  \caption{White light interferometry measurement of the topography of the Corona Heights fault plane at the millimeter scale (\cite{fault_database, doi_fault_database}). This surface, called $S_{mm}$, has a radius of $1.5$\,mm and is defined on a $2 \,\upmu$m grid with a $0.025\,\upmu$m elevation resolution. A parametric angle $\theta$ is defined from the main groove orientation.}
  \label{corona_millimetric}
\end{figure}
\begin{figure*}
  \includegraphics[width=0.8\linewidth]{PLAfault2.png}
  \caption{(Left): 3D file (STL file) obtained from the measured topography of the $S_{m}$ fault. The four 3D files generated for our friction experiments ($S_{m}$, $S_{mm}$ and their respective complementary) are available as Supplementary Material. (Right): Picture of a Polylactic acid 3D-print of the $S_{m}$ fault.}
  \label{PLAfault}
\end{figure*}
\\In seismic faults, such preferential orientations in their topography are observed at all scales (\cite{fault_imaging, fault_imaging2, rock_surf_bandwidth}) and originate from several processes. At the molecular level, rock forming crystals may display some frictional anisotropy. It is notably the case for antigorite, a mineral abundant in the Earth's upper mantle (\cite{friction_anisotropy_serp}). At the mesoscopic scale, the shear strength of foliated rocks is known to be anisotropic, due to the oriented planes in their constitutive mineralogy (\cite{strength_aniso1, strength_aniso2}). Fault zones in sedimentary basins are initiated by early fractures that often propagate in layered sediments. It can result in an anisotropic ramp-flat morphology of these fracture surfaces (\cite{rock_anisotropy}). For more mature faults having accumulated enough displacement, and above a given length scale (\cite{minimum_grooves}), the topography of the fault planes is also marked by slip induced wear, with striations and grooves of various wavelengths and amplitudes oriented along the main direction of slip (\cite{micro_grooves, macro_grooves}). If such morphological anisotropy of fault surfaces is well-known, its effect on the anisotropy of the frictional properties remains to be characterised.\\
Such a characterisation of frictional anisotropy could also be of interest for other types of rock contacts than strictly seismic faults, in particular for shallow rock joints and fractures, whose three-dimensional geometry is key in geotechnical engineering and for the structural stability of many man-made constructions (\cite{joint4, joint3, joint2, joint1}).\\
Here, we study how the morphology of faults controls the static coefficient of friction and the anisotropy of friction with regards to the main stress orientation during slip. To reach this goal, we produce 3D-prints of actual faults surfaces whose topography was measured in the field (\cite{fault_database}). We perform friction experiments with plaster of Paris casts of these 3D-printed faults. Results show that the coefficient of static friction along faults is highly anisotropic, a property that should henceforward be considered in numerical models of slip on seismic faults. We also show that this anisotropy is stress dependent, and should decrease with depth (e.g.,\,\cite{Byerlee1978}).

\section{3D printing and plaster casting of fault planes}

The actual morphology of natural faults can be difficult to assess, even if their long wavelength structures can be inferred by surface or subsurface imaging techniques (\cite{macro_grooves, imaging_sat, imaging_mag}). Yet, some fault planes are accessible to direct, high resolution, measurements, notably as they were exhumed by erosion and tectonic processes. For this study, we have used a series of digital fault surfaces. These fault roughness data were acquired with Light Detection And Ranging (LiDAR), laboratory laser profilometry, or white light interferometry techniques (\cite{fault_imaging2}). These data are available on an online public database (\cite{fault_database}), and in a repository with a doi number (\cite{doi_fault_database}). Should the reader hold some similar data, these authors welcome additions to this database. We have specifically selected fault roughness measurements performed on the Corona Heights fault (\cite{fault_imaging_3}) that outcrops near the Peixotto playground in San Francisco, California. These data cover surface areas with spatial scales in the range of millimeters to meters. Figures \ref{corona_metric} and \ref{corona_millimetric} show the fault surface at two spatial scales, one surface at the meter scale, defined on a $5\,$mm$\times5\,$mm grid, and one surface at the millimeter scale, defined on a $2\,\upmu$m$\times2\,\upmu$m grid. We will further on refer to these two surfaces as respectively $S_{m}$ and $S_{mm}$. Already, one can notice some preferential orientations in these topographies, and that the amplitude of fault roughness is, relatively to their size, somewhat larger at smaller scales\,($S_{mm}$) than at larger scales\,($S_{m}$) (\cite{fault_imaging_4}).
\begin{figure}[b]
  \includegraphics[width=1\linewidth]{plasterfaults2.png}
  \caption{Plaster faults made using the 3D-printed moulds (i.e., an example of which is shown in Fig.\,\ref{PlasFaults}). Left: $S_{m}$ (top) and its complementary surface (bottom). Right: $S_{mm}$ (top) and its complementary surface (bottom). The samples have a diameter of $40$\,mm. One can appreciate that the fault at the millimeter scale ($S_{mm}$) shows a higher roughness aspect than at the meter scale ($S_{m}$).}
  \label{PlasFaults}
\end{figure}
\begin{figure*}
  \includegraphics[width=0.9\linewidth]{setup2.PNG}
  \caption{Schematic (left) and picture (right) of the shear apparatus. It contains: the complementary plaster-casted fault surfaces (1 and 2) installed between two clamps; the compression spring (3) necessary to obtain a null normal loading when required (not mounted on the picture); and two horizontal and vertical sliders (4) used to force the motion in the direction of interest while allowing for vertical displacement. The shear force $F_T$ is applied on the top fault wall in the $\boldsymbol{e_\theta}$ direction, while the bottom surface is kept fixed. It is measured by a Sauter\textsuperscript{\textregistered} FH500N force gauge (\cite{Sauter}). The normal force $F_N$ is applied by a dead weight that acts oppositely to the $\boldsymbol{e_Z}$ direction) on the top surface. For a visualisation purpose, the two plaster blocks (1 and 2) are offset vertically in the schematic, whereas in the experiments these two surfaces matched, with a quasi-null aperture.}
  \label{setup}
\end{figure*}
\\For our tests, we chose to limit these anisotropic surfaces to a circular sample geometry. We also applied a mild running-window median filter to smooth out spikes in the measured surfaces that could be associated to measurement noise. The window length of the filter was $10$ space steps, accounting for $5\,$cm for $S_{m}$ and for $20\,\upmu$m for $S_{mm}$. In order to run the friction experiments, we generated some opposing surfaces to the ones presented in Figs. \ref{corona_metric} and \ref{corona_millimetric}. These opposing surfaces could not be measured, as the actual fault walls that were facing $S_{m}$ and $S_{mm}$ are now eroded. To reconstruct them, we have applied the following transformation to the 3D coordinates $(X, Y, Z)$ of $S_{m}$ and $S_{mm}$:
\begin{equation}
\begin{split}
    &X' \sim X\\
    &Y' \sim -Y\\
    &Z' \sim -Z,
\end{split}
\end{equation}
where $X'$, $Y'$ and $Z'$ are the coordinates of the generated opposing surfaces and $(X, Y)$ give the map location of a given surface point of elevation $Z$, as represented in Fig.\,\ref{corona_metric}. We have thus assumed that the missing fault walls are complementary to the measured ones, so that, when pressed together before the friction tests, they form a bulk with negligible aperture between the two blocks. Such assumption for natural faults would only be partly verified. When having accumulated enough slip, a granular layer of gouge material may there have formed, and the two opposing sliding surfaces may not always perfectly match. However, our assumption is relevant for the youngest faults with a small amount of slip. We have also assumed that erosion did not significantly alter the fault plane, such that the measured topography is representative of the one of an actual buried fault. For the Corona Heights fault, this assumption is valid because the fault offsets silica-rich chert rocks with a high resistance to weathering.
\\After having obtained the surfaces, we isotropically (i.e., with the same factor in all directions) down- or up-scaled $S_{m}$ and $S_{mm}$ to fit a standard $4$\,cm diameter disk that matches the clamp size of our shear deformation apparatus. We also re-gridded the surfaces to match the lateral resolution of our 3D-printer (Ultimaker\textsuperscript{2} Extended+\,(\cite{Ultimaker})) that has a nozzle size of $250\,\upmu$m. The four surfaces (two fault surfaces and two opposing surfaces) were then 3D-printed into polylactic acid (PLA) material, as shown in Fig.\,\ref{PLAfault}.
It should be noted that, even when designed to be flat, printed objects can present a natural roughness (\cite{3dprint_topo}), at a scale however smaller that the grooves observed on the printed faults. These intrinsic imperfections shall be comparable to $60\,\upmu$m, the elementary thickness of the PLA layers deposited by our 3D printer. In comparison, the 2D standard deviation of the elevation in our printed objects topography are $0.66\,$mm for $S_{m}$ and $1.7\,$mm for $S_{mm}$. The maximal elevation of these objects are, respectively, $3.7$\,mm and $8.3$\,mm. We thus consider that the small scale roughness ($\sim 60\,\upmu$m) from the printer's limit in resolution has a second order effect on the frictional properties of the surfaces.\\
Although we could have performed the friction experiments with the plastic pieces produced with the 3D printer, we have rather produced samples of plaster of Paris (gypsum) blocks moulded from the plastic faults. Plaster is known to be a reasonable model of porous brittle materials (\cite{plaster_modeling}), and the main goal of these casts was to work with a rock-like material, notably because plaster may wear and deform differently than plastic under shear. The fragile nature of plaster, and the potential friction-induced wear that the plaster was subjected to in our experiments, made us use new casts for each experimental realisation. The casts were generated with the following protocol: five volumes of water and eight volumes of powder of plaster of Paris were mixed and poured over the plastic moulds, then let to dry during one and half hour. The moulds, an example of which is shown in Fig.\,\ref{PLAfault}, were sprayed before each cast with a thin layer of silicon grease to avoid some of the fine plaster details to stick to the plastic during the mould release. The last step in the casts preparation was to dry them in an oven at a temperature of $40^\circ$C for one hour. As a result, we produced fault planes in blocks made of plaster of Paris, as shown in Fig.\,\ref{PlasFaults}.

\section{Experimental set-up and experimental conditions}

The shear apparatus used to perform the friction tests is shown in Fig.\,\ref{setup}. The two complementary surfaces are pressed together and mounted one on top of the other between the clamps of the shear apparatus. A normal force $F_N$ is applied on the top surface by using adjustable weights. In addition, a spring system of stiffness $625$\,N\,m\textsuperscript{-1} allows, if desired, to compensate for the machine empty weight of $13.7$\,N (i.e., the normal weight transmitted to the friction surfaces by the machine top clamp and structure and the top cast when no extra mass is used). A tangential driving shear force $F_T$ is then applied to the top fault wall in a given direction of the $(X, Y)$ plane. The amplitude of the force is measured by a Sauter\textsuperscript{\textregistered} force gauge (\cite{Sauter}).
The shear orientation is defined by the angle $\theta \in [0^\circ\,\,360^\circ]$ from the orientation of main grooves on the fault surface, as defined in Figs. \ref{corona_metric} and \ref{corona_millimetric}. A horizontal mechanical slider makes sure that the friction is evaluated in the direction of interest only, and a vertical slider allows upward or downward displacement of the top surface. While these sliders would ideally be perfectly lubricated, we have estimated their frictional resistance at the sliding velocity of our experiments, $F_c = 4\,$N$\pm0.5\,$N by performing a friction test with no fault installed in the machine (that is, with only air between the two clamps represented in Fig.\,\ref{setup}).
\begin{figure}[b]
  \includegraphics[width=1\linewidth]{force_time.PNG}
  \caption{Typical tangential force versus time for a given slide angle ($\theta=30^\text{°}$) and a given normal force ($F_N=33.3$\,N) applied to the $S_m$ fault. From the local maximal value of the tangential force, at the onset of slip, a friction coefficient $\mu_s(33.3\,\text{N},\,30^\text{°})=2.0$ can here be calculated using Eq.\,(\ref{coulomb}).}
  \label{courbe}
\end{figure}
\\
The target speed of the test bench speed (that is, the demanded slip velocity) was fixed to a constant and equal to $1.3$\,mm\,s\textsuperscript{-1}. Of course, such a velocity may be orders of magnitude above that of typical tectonic solicitations (for instance, an ultra-fast oceanic ridge may reach an opening rate of $10$ to $20$\,cm\,yr\textsuperscript{-1} (\cite{Fast_ridge})). Here, we define the static friction in the experiments as the peak shear stress reached before sliding occurs, divided by the normal stress.\\
We characterise the anisotropy of this laboratory static friction coefficient for both the $S_{m}$ and $S_{mm}$ surfaces by performing a series of experiments where we vary the angle of loading with respect to the grooves. We ran friction tests every $30^\circ$ on both plaster faults. At each angle, the experiment was repeated at least three times in order to ensure that results are reproducible (with new casts each time, to avoid any wear related deviation). The standard deviation computed on these multiple measurements (typically $5$ to $10$\,N) was used to compute the error bars on our characteristics coefficients of static friction. $S_m$ was sheared under a normal stress $\sigma_N=10.9$\,kPa, while the tests performed on the rougher surface $S_{mm}$ were performed under $\sigma_N=26.5$\,kPa. A total of 76 experiments were performed, 37 using $S_{m}$ and 39 using $S_{mm}$.
At the onset of slip, the laboratory static friction coefficient is defined using a standard Coulomb's law (\cite{Bowden1951}):
\begin{equation}
    \mu_s(F_N, \theta) = \cfrac{\max\left[F_T(F_N, \theta)\right]-F_c}{F_N},
    \label{coulomb}
\end{equation}
where $\mu_s$ is the coefficient of static friction and $\max\left[F_T\right]$ is the (local) maximum tangential force applied at the onset of slip. In the following, we will also consider the mean driving and normal stresses, denoted $\sigma_T=F_T/(\pi r^2)$ and $\sigma_N=F_N/(\pi r^2)$, where $r=2\,$cm is the radius of the cast. Figure \ref{courbe} shows a typical measurement of a friction test, from which $\max\left[F_T\right]$, and hence $\mu_s$, are calculated.

\section{Results}

\subsection{Friction anisotropy}

The results are presented in Fig.\,\ref{results}. The derived coefficient of frictions are larger than $1$, which does not come as a surprise due to the non negligible roughness of our fault samples. Indeed, a large part of the contact area is bound to be perpendicular to the demanded slip (in any direction), inducing a strong resistance to motion.\\
On both fault samples, one can observe the strong anisotropy of the coefficient of static friction, with the maximum value of $\sigma_S$ being about four times larger than its minimum for $S_m$ and about three times larger for $S_{mm}$. In most experiments, the minimal friction is obtained along the main groove orientation (i.e., at $\theta=0^\circ$ or $\theta=180^\circ$), and the resistance to shear is larger perpendicularly to this orientation. The maximum of value is however never obtained exactly at $\theta=90^\circ$ or $\theta=270^\circ$ but rather along a neighbouring direction. Local maxima are indeed obtained for $\theta=60^\circ$ or $\theta=240^\circ$ when shearing $S_{m}$ and for $\theta=120^\circ$ or $\theta=300^\circ$ when shearing $S_{mm}$.

\begin{figure}
  \includegraphics[width=1\linewidth]{sm_smm_mu.PNG}
  \caption{Coefficient of static friction of the $S_m$ fault (top) and $S_{mm}$ fault (bottom) as a function of the sliding direction $\theta$. Each transverse error bar is computed using at least three experimental realisations (with new plaster casts on each occurrence). Bar lengths display twice the standard deviation of the obtained coefficient. The tests performed with sample $S_m$ were done under a normal stress $\sigma_N=10.9$\,kPa, while those performed with the sample $S_{mm}$ (larger roughness) were done under $\sigma_N=26.5$\,kPa. The arrows indicate the orientation of the main grooves (see Figs.\,\ref{corona_metric} and.\,\ref{corona_millimetric}).}
  \label{results}
\end{figure}

\subsection{Damage and stress dependence}

\begin{figure}
  \includegraphics[width=1\linewidth]{damages.PNG}
  \caption{Various types of damage observed on the fault surfaces after the sliding tests. The dashed arrows show the slip directions of the complementary plaster casts during the tests. (a): $S_m$ sheared at $\theta=330^\circ$ and $\sigma_N=26.5$\,kPa. The indicated broken topographic high, on top of the fault, comes from the complementary surface. (b): Complementary of $S_m$ of a different experiment, sheared at $\theta=90^\circ$ and under a higher normal stress $\sigma_N=99.0$\,kPa. The visible damage is there less localised. Scale: the samples diameter is $40\,$mm.}
  \label{damages}
\end{figure}
\begin{figure}
  \includegraphics[width=1\linewidth]{profile_60240.PNG}
  \caption{Elevation profile along the diameter of $S_m$ in the $60^\circ-240^\circ$ direction. The distances are reported in real word units and the axes are orthonormal. The arrows indicate the motion of the complementary surface for tests along $\theta=60^\circ$ and along $\theta=240^\circ$. The chip is only strongly interlocked in the latter orientation.}
  \label{profile}
\end{figure}
\begin{figure}
  \includegraphics[width=1\linewidth]{stress_dep.PNG}
  \caption{Coefficient of static friction $\mu_S$ as a function of the applied normal stress $\sigma_N$ for the fault $S_m$ sheared along the direction $\theta=90^\circ$ (squares). The length of the error bars is twice the standard deviation of $\mu_S$ obtained for three different tests. For reference, the straight dashed line indicates $\mu_S\propto {\sigma_S}^{-0.45}$. The static coefficient of friction measured on contacting flat plaster surfaces is shown (crosses) for comparison.}
  \label{stress_dep}
\end{figure}
Most of our experiments were destructive, with visible wear on the plaster samples after the shearing tests. This wear was the main reason calling for the production of new plaster casts for each experimental realisation, as we verified that repeating a same experiment with a previously used cast led to a significant (and here unwanted) drop in friction. The observed damage consists either in the formation of plaster powder (gouge) or in the rupture of topographic highs of the fault surfaces. Part of it might have initiated at the onset of the fault displacement (and hence be related to the static friction), while some of it has rather been induced by the subsequent sliding. Figure\,\ref{damages} shows some examples of these damage types.\\
Wear of seismic faults has been studied (e.g.,\,\cite{fault_wear, fault_damage_zones}) to, in particular, better understand the energy budget of the deformation, but also because this process may lubricate faults during slip (e.g.,\,\cite{fault_lubri,FluidPresMelt}) or modify the fault permeability to fluid flow (e.g.,\,\cite{gouge_flow1,gouge_flow2}). The present study focuses on the measurement of the coefficient of static friction and on its anisotropy, but we suggest that our 3D-print-based set-up could also enable the quantitative characterisation of damage during sliding along analogue fault surfaces.\\
We here keep to a qualitative assessment of which parts of the surfaces were mainly worn during each experiment. It seems that most of the shear resistance of the Corona Heights fault, at the millimeter scale ($S_{mm}$), arises from its grooves. By contrast, the friction of $S_{m}$ (representing a metric scale) is dominated by one chip in its field-scanned morphology (i.e., the main topographic low in Fig.\,\ref{corona_metric}). This chip being on the edge of the 3D-printed surface, but not on the edge of the real word fault, a finite size effect (inducing an artificial asymmetry of the fault wall) is certainly at play in the results reported in Fig.\,\ref{results} (top), notably explaining the strong asymmetry in $\mu_S$ for the opposite directions $\theta=60^\circ$ and $\theta=240^\circ$. Indeed, due to the finite size of our tested pieces, the interlocking of the chip is relatively free to unlock laterally along the $\theta=60^\circ$ direction, but is strongly locked by the surrounding plaster along the $\theta=240^\circ$ direction (see Fig.\,\ref{profile}).\\
Because the overall friction is likely to be affected by the surface wear, and because this wear is likely stress dependent, we have performed some friction tests on $S_m$ under various loads $\sigma_N$, at a given angle $\theta=90^\circ$. The results are shown in Fig.\,\ref{stress_dep}. At the highest tested stresses, $\mu_S$ seems, to an extent stress independent, with its mean variations lying within the measured error bars. While this result is compatible with the classical Coulomb theory (e.g.,\,\cite{Bowden1951}), one can observe, over a wider range of normal stresses, a consistent decrease in the friction coefficient with a higher normal stress. It could, in part, emanate from some limitations in our experimental set-up. For instance, at lower $\sigma_N$, the internal friction of the device ($F_c$) accounts for a more significant portion of the total measured tangential force, and the friction characterisation could thus be less accurate. We have however run similar tests on flat plaster surfaces, showing no significant variations in $\mu_S$ over the same stress range (see Fig.\,\ref{stress_dep}). The drop in friction coefficient with $\sigma_N$ is then likely related to the increased damage under higher normal loading (see Fig.\,\ref{damages}), reducing the overall shear resistance as asperities are easier to break. It might also result from the change in effective contact area with a higher load, with, proportionally, more pressure being borne by the surfaces parallel to the demanded motion, causing proportionally less of a resistance to slip. Such a drop of $\mu_S$ with the normal stress has already been reported for tilted contacts (\cite{Oded_Fineberg}), as the static coefficient of friction is not an absolute material constant (\cite{Oded_Fineberg, fractal_friction}).\\
We have then assessed the effect of the normal stress $\sigma_N$ on the anisotropy of the static friction coefficient. We performed frictional tests on the four poles of $S_m$ ($\theta=0^\circ,\,90^\circ,\,180^\circ \text{ and } 270^\circ$) at a higher load ($\sigma_N=171$\,kPa) than the load used before (i.e., $\sigma_N=10.9$\,kPa, as reported in Fig.\,\ref{results}). The newly measured coefficients of static friction are shown in Fig.\,\ref{aniso_stress}. One can notice the reduced friction anisotropy at high $\sigma_N$. The ratio between the maximum and the minimum value of $\mu_S$ indeed drops from $3.1$ at $\sigma_N=10.9$\,kPa to $1.5$ at $\sigma_N=171$\,kPa. This result suggests that the frictional anisotropy of faults is smaller at depth. This concept is naturally linked with the seminal \citet{Byerlee1978}'s law, stating that the roughness of fault planes (as well as the type of their constitutive rocks) has less effect on their maximum static frictional properties at larger depths.
\begin{figure}
  \includegraphics[width=1\linewidth]{aniso_stress.PNG}
  \caption{Coefficient of static friction $\mu_S$ of the $S_m$ fault along the cardinal directions $\theta=0^\circ,\,90^\circ,\,180^\circ \text{ and } 270^\circ$, for two normal stresses: $\sigma_N=10.9\,$kPa (outer plain line) and $\sigma_N=171\,$kPa (inner dashed line). The bars length is twice the standard deviation of the measured coefficient, for three different experimental realisations, at each angle and stress. Friction data at the two different normal loads do not overlap. One can notice the consistent reduction in the friction coefficient with a higher load and the reduced friction anisotropy at high $\sigma_N$.}
  \label{aniso_stress}
\end{figure}

\section{Discussion and conclusion}

Here, we show how the multiscale anisotropy of fault plane topography leads to an anisotropy in the frictional properties. Results confirm that seismic faults are prone to slide along some preferential orientations. The orientation that is the most likely is the one that faults have previously slid along, and which has shaped some guiding grooves in their morphology. Yet, displacements following other orientations are possible. Predicting the rupture direction of the next earthquake on a fault is thus not only dependent on assessing the main regional stress. The question should rather be along which orientation a rupture criterion (e.g., \cite{Bowden1951}) will first be exceeded. Such a subtlety might be of little importance for mature faults for which the stress principle orientations have not changed with time, because, in this case, the main stress is likely to act along the lowest coefficient of friction anyway. Yet, it could be paramount for faults under a changing geological load, where this alignment is not verified, or for immature faults, where the slip could be mainly governed by the anisotropy of early surfaces (i.e., where the slip does not coincide with the stress principal orientation, but is non-associated). Examples of slickensides (i.e., fault planes) commonly exist with several overlapping striations orientations (e.g.,\,\cite{multiple_orientations}), with rake and striations oblique to the actual orientation of the fault plane\,(e.g.,\,\cite{oblique_striations}). These observations indicate that the original slip direction (if assumed to initiate following Andersonian criterion (\cite{Anderson_faulting})) does not completely determine the direction of the next episode. Earthquakes occurring along abnormal directions (i.e., not in agreement with the local stress state) have been observed (e.g.,\,\cite{Balochistan_earthquake, tensile_Torishima_EarthQ, Anderson_faulting}), and their understanding might be eased by accounting for the possible frictional anisotropy of their surfaces (\cite{oblique_slip, stress_inversion_aniso}).\\
Note that frictional anisotropy should not only be considered at the full fault scale, but this property may spatially vary along the fault walls. Analytical solutions demonstrate that the stress around a fault is perturbed by its roughness, and a local slip can occur much before the entire fault is under yielding conditions in a given direction (\cite{interlocked_rough_faults}).\\
While we have here only measured the static coefficient of friction, we suggest that similar studies could be performed to characterise the coefficient of dynamic friction (i.e., by analysing the evolution of the resistance to motion, after the plaster faults start moving, as a function of the sliding orientation). Hence, not only the initial slip direction of an earthquake could be impacted by frictional anisotropy, but the complete slip trajectory (\cite{friction_anisotropy_2}). Changes in the slip direction within single earthquake rupture events are indeed sometimes observed, notably from bent grooves on post-mortem fault walls (\cite{variable_rake}).\\
We have, additionally, measured how the anisotropy in friction becomes less significant when the normal stress acting on a fault increases (i.e., with the fault depth), in general agreement with Byerlee's law (\cite{Byerlee1978}).
Such an effect likely derives from the stress related changes in rupture rheology and in damage type. The transition from a highly anisotropic to a relatively isotropic regime should typically occur when local stresses on the fault reach the yield strength of the material, $\sigma_y\sim5$\,MPa in the case of plaster (\cite{plaster_modeling}). This is about two orders of magnitude above the transition $\sigma_N\sim100\,$kPa at which we observed a strong reduction in anisotropy (see Fig.\,\ref{aniso_stress}), but our computed $\sigma_N$ is an average value which does not account for the potential strong stress concentrations at play in our faults.
Considering that the strength $\sigma_y$ of rocks (e.g.,\,\cite{rock_strength}) is about two orders of magnitude (100 times) larger than that of plaster, fault frictional anisotropy could thus be only at play at pressures less than about $\sigma_N\sim100\times100\,\text{kPa}=10\,$MPa. This would correspond to the shallowest faults, at depths less than $\sigma_N/(g\rho)\sim500\,$m, where $\rho$ is the volumetric mass of rocks ($\sim2000\,$kg\,m\textsuperscript{-3}) and $g$ the gravity acceleration. Some care should however be taken when deriving such a conclusion by analysing resized samples as ours ($40\,$mm diameter samples representing meter or millimeter topographies), as the way matter breaks is length-scale dependent\,(e.g.,\,\cite{fault_imaging_4}).
Note also that other fault geometries than the one we have here studied may induce lesser stress concentration, so that a significant damage only occurs at a mean stress level directly comparable to the yield stress $\sigma_y$ of rocks. Thus, the possibility of frictional anisotropy should not be overlooked when studying fault buried up to $\sigma_y/(g\rho)\sim50\,$km. The heterogeneity of fault planes, and thus the anisotropy in this heterogeneity, may also still play a role under high stress, as roughness does not only encourage local yield, but also helps to suppress large slip events on moving faults (\cite{roughness2}).\\
Additionally to the assessment of the stability of (at least) shallow seismic faults, the characterisation of the frictional anisotropy of rock surfaces may be of importance in geotechnical engineering, for instance, for the stability of tunnels and foundations. There, the intrinsic strength anisotropy of foliated rocks is well studied (\cite{strength_aniso1, strength_aniso2}). Our work shows how one can also characterise the mechanical anisotropy of rough rock contacts, for instance, along joints (\cite{joint4, joint3, joint2, joint1}) and fractures (\cite{fracture1, fracture2}) between or inside rock formations.\\
A main point of this manuscript is, finally, to illustrate how the 3D-printing technology can help with new experimental designs in Earth Sciences, and this technology is getting a growing attention from the community (\cite{3dprint1, 3dprint2, 3dprint3, 3dprint_poro,3dprint_fault}), including the study of the frictional properties of 3D-printed fault analogues (\cite{3dprint_fault}). A direct continuation of the present work, for instance, could be to 3D-print and to test some faults surfaces beforehand filtered with various band-pass filters, in order to understand how the various wavelengths of the topography contribute to the global static friction coefficient, to the dynamical friction coefficient and to analyse the spatial distribution of the fault wear produced under various stresses and amounts of slip.\\

\subsection*{Acknowledgement}

\noindent
\textbf{Author contributions:} RT proposed the guidelines of this work, AS built the test machine and, with TVD, printed the faults and performed the friction experiments. FR advised on the mechanics of seismic faults. TVD wrote the first version of this manuscript and all authors contributed to the writing of its final version. We are grateful for the early experimental explorations performed by Marine-Sophie Jacob, Céline Fliedner, Aldo Mellado Aguilar, Gaëtan Leca and Laifa Rahmi, all students from the EOST/IPGS faculty at the University of Strasbourg. We also thank Amir Sagy from the Geological Survey of Israel for fruitful discussions, and acknowledge the support of the University of Strasbourg, of the IRP France-Norway D-FFRACT, and of SFF Porelab (project number 262644 of the Research Council of Norway). We thank the Strasbourg AV.Lab association for the use of their 3D printer.
\textbf{\\Author declarations:} The authors declare no competing interest. A funding support from the University of Strasbourg is acknowledged. Readers are welcome to comment and correspondence should be addressed to vincentdospitalt@unistra.fr or renaud.toussaint@unistra.fr.

\newpage
\bibliographystyle{plainnat_for_frontiers.bst}