Fix rendering docs to PDF (#735)

I've been fixing some problems with rendering the docs to PDF. You can
try it yourself by running:

`quarto render docs --to pdf`

I've hit a roadblock with the following error which first appears when
attempting to render `modflow-demo.qmd`:
`Unable to load picture or PDF file
'https://user-images.githubusercontent.com
/13662783/187665858-d01fd60f-f3c2-4662-af82-cf8acfbe169b.PNG'.`

It seems that [_.. pdfTeX and XeTeX do not include the necessary code to
grab an image from an arbitrary location.
..._](https://tex.stackexchange.com/questions/5433/can-i-use-an-image-located-on-the-web-in-a-latex-document)

Maybe a workaround is possible with pixi that downloads the necessary
files to a local folder?

docs/core/allocation.qmd

@@ -172,7 +172,7 @@ $$
 $$
     \hat{C}_S^p \leftarrow \hat{C}_S^p - F;
 $$
-5. Repeat steps 2-4 for the remaining priorities up to $p_\max$.
+5. Repeat steps 2-4 for the remaining priorities up to $p_{\max}$.
 
 :::{.callout-note}
 In the future there will be another solve before the priority 1 solve, taking the demand/supply from basins into account.

docs/core/equations.qmd

@@ -94,18 +94,16 @@ The presence of division by the basin area means that areas of size zero are not
 ## The reduction factor {#sec-reduction_factor}
 At several points in the equations below a *reduction factor* is used. This is a term that makes certain transitions more smooth, for instance when a pump stops providing water when its source basin dries up. The reduction factor is given by
 
-$$
+\begin{align}
     \phi(x; p) =
-    \begin{align}
-        \begin{cases}
-        0 &\text{if}\quad x < 0 \\
-            -2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\
-            1 &\text{if}\quad x > p
-        \end{cases}
-    \end{align},
-$$
+    \begin{cases}
+    0 &\text{if}\quad x < 0 \\
+        -2 \left(\frac{x}{p}\right)^3 + 3\left(\frac{x}{p}\right)^2 &\text{if}\quad 0 \le x \le p \\
+        1 &\text{if}\quad x > p
+    \end{cases}
+\end{align}
 
-where $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.
+Here $p > 0$ is the threshold value which determines the interval $[0,p]$ of the smooth transition between $0$ and $1$, see the plot below.
 
 ```{python}
 # | code-fold: true
@@ -463,23 +461,21 @@ $$
     Q_\text{PID}(t) = K_p e(t) + K_i\int_{t_0}^t e(\tau)\text{d}\tau + K_d \frac{\text{d}e}{\text{d}t},
 $$ {#eq-PIDflow}
 
-for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_\min, Q_\max$, and so
+for given constant parameters $K_p,K_i,K_d$. The pump or outlet can have associated minimum and maximum flow rates $Q_{\min}, Q_{\max}$, and so
 $$
-Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_\min, Q_\max).
+Q_\text{pump/outlet} = \text{clip}(\Phi Q_\text{PID}; Q_{\min}, Q_{\max}).
 $$
 
 Here $u_\text{us}$ is the storage of the basin upstream of the pump or outlet, $\Phi$ is the product of [reduction factors](equations.qmd#sec-reduction_factor) associated with the [pump](equations.qmd#sec-pump) or [outlet](equations.qmd#sec-outlet) and
 
-$$
-\text{clip}(Q; Q_\min, Q_\max) =
 \begin{align}
+    \text{clip}(Q; Q_{\min}, Q_{\max}) =
     \begin{cases}
-        Q_\min & \text{if} & Q < Q_\min \\
-        Q & \text{if} & Q_\min \leq Q \leq Q_\max \\
-        Q_\max & \text{if} & Q > Q_\max
+        Q_{\min} & \text{if} \quad Q < Q_{\min} \\
+        Q & \text{if} \quad Q_{\min} \leq Q \leq Q_{\max} \\
+        Q_{\max} & \text{if} \quad Q > Q_{\max}
     \end{cases}.
 \end{align}
-$$
 
 For the integral term we denote
 $$
@@ -511,7 +507,7 @@ that is, $\hat{f}_\text{PID}$ is the right hand side of the ODE for the controll
 
 Using this, solving @eq-PIDflow for $Q_\text{PID}$ yields
 $$
-Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_\min,Q_\max\right),
+Q_\text{pump/outlet} = \text{clip}\left(\phi(u_\text{us})\frac{K_pe + K_iI + K_d \left(\frac{\text{d}\text{SP}}{\text{d}t}-\frac{\hat{f}_\text{PID}}{A(u_\text{PID})}\right)}{1\pm\phi(u_\text{us})\frac{K_d}{A(u_\text{PID})}};Q_{\min},Q_{\max}\right),
 $$
 where the clipping is again done last. Note that to compute this, $\hat{f}_\text{PID}$ has to be known first, meaning that the PID controlled pump/outlet flow rate has to be computed after all other contributions to the PID controlled basin's storage are known.
 

docs/core/numerics.qmd

@@ -49,25 +49,24 @@ $$
 \mathbf{g}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1})
 $$
 by iterating $\mathbf{g}$ on an initial guess of $\mathbf{w}_{n+1}$;
-- Newton-Raphson iterations: approximate $\mathbf{w}_{n+1}$ as a root of the function
+- Newton iterations: approximate $\mathbf{w}_{n+1}$ as a root of the function
 $$
 \mathbf{h}(\mathbf{x}) = \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{x},t_{n+1}) - \mathbf{x},
 $$
 by iteratively finding the root of its linearized form:
-$$
+
 \begin{align}
 \mathbf{0} =& \mathbf{h}(\mathbf{w}_{n+1}^k) + \mathbf{J}(\mathbf{h})(\mathbf{w}_{n+1}^k)(\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k) \\
 =& \mathbf{w}_n + (t_{n+1}-t_n)\mathbf{f}(\mathbf{w}_{n+1}^k,t_{n+1}) - \mathbf{w}_{n+1}^k \\ +&\left[(t_{n+1}-t_n)\mathbf{J}(\mathbf{f})(\mathbf{w}_{n+1}^k)-\mathbf{I}\right](\mathbf{w}_{n+1}^{k+1}-\mathbf{w}_{n+1}^k).
 \end{align}
-$$ {#eq-newtoniter}
 Note that this thus requires an evaluation of the Jacobian of $\mathbf{f}$ and solving a linear system per iteration.
 
 # The advantage of implicit methods
 
 The implicit method @eq-eulerbackward is a coupled system of equations for $\mathbf{w}_{n+1}$, while the explicit @eq-eulerforward is fully decoupled. This means in general that in a time integration step with an implicit method the basins communicate information with eachother, while in an explicit method this does not happen. A consequence of this is that local events (e.g. a pump turns on) propagate slowly trough the model using an explicit method but quickly using an implicit method, making implicit methods more stable.
 
 # Jacobian computations
-The iterations @eq-newtoniter require an evaluation of the Jacobian of $\mathbf{f}$. The Jacobian of `water_balance!` is discussed [here](equations.qmd#the-jacobian).
+The Newton iterations above require an evaluation of the Jacobian of $\mathbf{f}$. The Jacobian of `water_balance!` is discussed [here](equations.qmd#the-jacobian).
 
 There are several ways to compute the Jacobian:
 
@@ -77,7 +76,7 @@ There are several ways to compute the Jacobian:
 
 # Continuity considerations
 
-The convergence of the Newton-Raphson method can be [proven](https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method) given certain properties of $\mathbf{f}$ around the initial guess and root to find. An important aspect is the smoothness of $\mathbf{f}$. The basin profiles and $Q(h)$ relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.
+The convergence of the Newton method can be [proven](https://en.wikipedia.org/wiki/Newton%27s_method#Proof_of_quadratic_convergence_for_Newton's_iterative_method) given certain properties of $\mathbf{f}$ around the initial guess and root to find. An important aspect is the smoothness of $\mathbf{f}$. The basin profiles and $Q(h)$ relations are given by interpolated data, and thus we have some control over the smoothness of these functions by the choice of interpolation method. This is discussed further below. the Manning resistance is not mentioned here since it is given by an analytical expression.
 
 Control mechanisms can change parameters of $\mathbf{f}$ discontinuously, leading to discontinuities of $\mathbf{f}$. This however does not yield problems for the time integration methods in `DifferentialEquations.jl`, since the [callback mechanisms](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/) used to change these parameters make the solver take these discontinuities into account.