07_BC3c.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example 7.\n",
    "### Sensitivity of hydraulic head at a point to <font color='blue'>Cauchy BC conductance</font> under steady state flow conditions\n",
    "\n",
    "## 0. Forward model\n",
    "<br>\n",
    "Governing equation:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad K\\ b\\ \\frac{d^2h}{dx^2}+R=0\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Boundary conditions:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad -K\\ b\\ \\frac{dh(x)}{dx}=0\\ ,\\quad\\quad\\quad\\quad\\quad\\quad\\;\\, x=0\\cr\n",
    "\\end{align}$\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad -K\\ b\\ \\frac{dh(x)}{dx}=C\\ [h^*-h(x)]\\ ,\\quad\\quad x=L\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Closed-form solution:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\cr\\end{align}$\n",
    "\n",
    "Spatial derivatives from differentiation:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import HTML, display\n",
    "def set_background(color):\n",
    "    script = (\n",
    "        \"var cell = this.closest('.code_cell');\"\n",
    "        \"var editor = cell.querySelector('.input_area');\"\n",
    "        \"editor.style.background='{}';\"\n",
    "        \"this.parentNode.removeChild(this)\").format(color)\n",
    "    display(HTML('<img src onerror=\"{}\">'.format(script)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "from warnings import filterwarnings\n",
    "filterwarnings(\"ignore\", category=DeprecationWarning)\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "K, R, b, L, BC3h, BC3c, ocol = 10., 1e-1/1000., 10., 10000., 1., 1., 5000\n",
    "X = np.arange(L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Direct sensitivity\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Perturbation sensitivity\n",
    "$\\begin{align}\\cr\n",
    "\\frac{\\partial h(x')}{\\partial C_{\\Gamma_3}}\\approx\\frac{h(x,C_{\\Gamma_3}+\\Delta C_{\\Gamma_3})-h(x,C_{\\Gamma_3})}{\\Delta C_{\\Gamma_3}}\\;\n",
    "\\cr\\cr\\end{align}$\n",
    "\n",
    "### 2a. Analytical\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Semi-analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAACdCAYAAACkc0HoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnD0lEQVR4nO3deVxVdf748dcHBBEXEHFFBDRXwBXTUkwFl9xySS211EbLaXImna9mNU1m06TpVKNpTtZPc6kJNRtNm0lFMms0UfkKhjtLuABioiIIwvv3B3S/rgQKHC68n4/Hfcg959zPed8jl/f9LOfzMSKCUkoppeyPg9UBKKWUUuruaBJXSiml7JQmcaWUUspOaRJXSiml7JQmcaWUUspOaRJXSiml7FQVqwMoLk9PT/H19bU6DKWUUqpM7Nu375yI1L3dPrtL4r6+vkRGRlodhlJKKVUmjDEJd9qnzelKKaWUndIkrpRSStkpu2tOL2l//OMfOXToEA4ODrZHixYtWLBgAQAzZ84kMTERBwcHHB0dcXBwoE2bNrzwwgsAvPnmm6Snp+Ps7IyzszNVq1alRYsWPPLIIwCsX7+e3Nxc2z5XV1e8vLy47777ADh37hzVqlWjWrVqODjodyqllFJFV+mTeEZGBhcuXCAvL8/2qFWrlm3/sWPHiI2NJTc317Y/JyfHtv+TTz7h2LFjXL161bZt6NChtiQ+ZcoUzp07d8M5x40bx6pVqwBo0qQJmZmZAFSrVo3q1aszefJk/vrXv5Kbm0vv3r2pUaMGtWrVws3NDTc3N0JDQ+nTpw85OTls2bIFNze3G/a7ubnh5ORUatdMKaVU+VDpk/jSpUsL3b9hw4ZC90dHRwMgIly7do3s7GyuX1Rmz549ZGZmkp2dTVZWFpmZmXh6etr2z58/n4yMDDIyMrhy5QoZGRkEBgYCkJOTgzGG5ORkjh49Snp6Ounp6bi4uNCnTx/S0tIYOnToLTHNnTuXF154gaSkJAYMGECdOnVsDw8PD0aMGEFQUBCXLl3i4MGDeHh42PZVqVLpfyWUUspu6F/sEmKMwcnJ6ZYacNOmTQt93e9+97s77nNxcSEiIuKW7Xl5eQB4eHgQGRlJeno6Fy9etCX5bt26AZCbm0uzZs1IS0sjNjaWtLQ00tLSaN26NUFBQURHR9O9e/cbynZzc2PVqlUMHjyYQ4cO8f7771O/fv0bHgEBAdSoUaMol0UppVQp0iRuh37pO3d2dqZTp053PM7Hx+eWlgQRsX0JaNOmDf/5z39syT0tLY1z587xy334iYmJrFmzhgsXLtxQxs6dOwkODmbdunXMmjWL+vXr07BhQ7y8vPDy8mLixInUrVuXjIwMHBwcqFatWsm9eaWUUjbG3tYTDwoKEr1PvGxdvXqVlJQUkpOTSU5Oplu3bri7uxMREcE//vEPkpOTOXPmDElJSVy+fJmTJ0/i5+fHW2+9xQsvvICHh4ctwXt5efG3v/0NNzc3Tp48SWZmJj4+PlqzV0qpOzDG7BORoNvu0ySuStLFixepUaMGDg4O/PDDD2zbto2kpCROnTrFqVOnOH36NPHx8Tg7O/P73/+eRYsWAVCnTh18fHzw8/MjLCwMBwcHYmJiuHbtGj4+Pri7u2OMsfjdKaVU2dMkrsql2NhYoqKiSEhIICEhgfj4eDIzM23jAIYNG8YXX3wBQM2aNfH19aVjx46sWLECgP3791O9enV8fX2pWrWqNW9CKaVKmSZxZZcOHTpEbGzsDUm+Vq1arF69GoBOnTqxf/9+jDE0adKE++67j5CQEF588UUAjh8/ToMGDbSpXill1wpL4jqwTZVb/v7++Pv733H/P/7xDw4fPszx48c5fvw4J06cIC4uzra/W7dupKSkUL9+fe677z5atGjBwIEDGTFiBJA/et/R0bHU34dSSpUWTeLKbgUFBREUdNsvp4gIixYt4sSJE7Yk/9VXX1G7dm1GjBhBZmYmtWvXxs/Pj1atWtke3bt3p1mzZmX8TpRS6u5oElcVkjGGUaNG3bI9NzcXgOzsbKZPn87hw4c5fPgwmzdvJicnh3feeYfnn3+ehIQExo8fT6tWrfD39ycgIICAgADq1r3taoBKKWUJTeKqUvml+dzNzY2//vWvtu3Xrl0jLi4ONzc3IH+UfU5ODmFhYfz888+249avX8/w4cM5fvw44eHhBAYG4u/vf8NUvUopVVY0iSsFVKlShebNm9ueBwYG8t133yEinD17lpiYGGJiYujcuTMA4eHhPPPMM7bjmzRpQkBAAEuWLMHHx4eLFy9StWpVHTWvlCpVOjpdqbuQl5dHQkIChw4dsiX46OhovvnmG9zd3fnzn//Mm2++SZs2bWjfvj3t27enQ4cOBAcH62A6pVSx6C1mSpWxXbt2sWXLFg4cOEBUVBRnz56lRo0apKen4+DgwPvvv09KSootwTdp0kQns1FK3ZbeYqZUGevevfsNi8ucPXuW+Ph427z34eHhrF+/3rbiXe3atRkyZIhtIpvk5GTq1auniV0pVShN4kqVgQYNGtCgQQPb87Vr15KRkUF0dDRRUVFERUXdsERtx44duXLlCkFBQXTu3JmgoCC6du1Ko0aNrAhfKVVOaXO6UuVMXl4eH330EZGRkezdu5fo6GiuXbvG1KlTWbhwITk5Obzxxhu2++Sv/3KglKp4tE9cKTuWlZXFwYMHcXNzo2XLlsTGxhIQEGBbUrZx48Z06dKF6dOn8+CDD1ocrVKqpGmfuFJ2zMXFhfvvv9/2vHXr1qSnpxMVFUVkZCQ//PADu3fv5vLly0B+f/usWbN44IEHeOCBB+jatSs+Pj7av65UBVSmSdwY4whEAqdEZJAxxgP4DPAF4oFRIvLznUtQSgHUqFHjlsFzv7SqiQiurq58+OGHLFy4EMjvk//+++/x8/Pj3LlzVK9enWrVqlkSu1Kq5DiU8fn+AMRe93wWsF1EmgPbC54rpe7CLzXtkJAQIiIiSE9PZ9++fSxevJj+/fvj7e0NwGuvvUatWrXo3Lkz06ZN4/PPPyclJcXK0JVSd6nM+sSNMY2Bj4E3gOkFNfEjQE8ROWOMaQhEiEjLwsrRPnGl7s2uXbv46quv+P7779m9ezdZWVk0atSIpKQkjDF8++23eHl54efnp03wSpUD5aVP/F1gJlDzum31ReQMQEEir1eG8ShVKV3fDJ+dnc3+/ftJSUmxJewnnniChIQEGjVqRHBwMN27dyc0NJRWrVpZGbZS6jbKpDndGDMISBGRfXf5+qeNMZHGmMjU1NQSjk6pysvZ2ZmuXbsyZMgQIL8/fdOmTSxZsoSHHnqI7777jqlTp/Luu+8C+avALViwgD179nDt2jULI1dKQRk1pxtj3gSeAK4BLkAt4HOgM9qcrlS5JSIkJCQgIvj5+REbG0ubNm0AqFmzJsHBwfTq1YuRI0fi4+NjcbRKVUyFNaeXSU1cRF4UkcYi4gs8BoSLyDhgIzC+4LDxwL/KIh6lVNEYY/D19cXPzw/Iv73t7Nmz/POf/2TMmDEcP36cGTNmcOjQIQBiYmL4+9//zsGDB233sSulSo/V94nPBcKMMb8BEoGRFsejlPoV9evXZ/To0YwePRqAU6dO4eHhAcDWrVuZPn06AJ6enjz00EP06tWLiRMn4urqalnMSlVUOmObUqpEJSYmsmPHDtsjNTWVCxcu4OzsTFhYGFevXqVPnz46XaxSRaTTriqlLCEiJCcn2xJ2SEgI4eHhAAQGBtK3b18GDRpEz549LYxSqfLN8j5xpVTlZIy5oca9detW9u3bx9y5c6lbty6LFi1i0aJFtv1Lly5l//792p+uVBFpTVwpZZkrV65w/vx5GjduzNmzZ2nYsCGQ358eGhpK3759GTBgAPXr17c4UqWsozVxpVS55OrqSuPGjYH8+d1Pnz7NypUrefjhh4mIiOCpp57iq6++AiA5OZmdO3eSk5NjZchKlStaE1dKlUsiQkxMDN7e3ri7u7N48WKee+453Nzc6Nu3Lw8//DD9+/e31d6Vqqh0YJtSyu5dvHiRbdu2sWXLFrZs2cKZM2dwcHDg3Llz1K5dm5SUFOrUqYOjo6PVoSpVosrL3OlKKXXXatWqxfDhwxk+fDgiwsGDB9m3bx+1a9cG8ud8j4yMtNXSBwwYgKenp8VRK1W6tCaulKoQNmzYwMaNG/nqq69ITk7GwcGBZ555hiVLllgdmlL3pERr4saY6kCWiOTec2RKKVVChg0bxrBhw8jLy+PAgQNs2rTJNl3spUuXePDBB+nXrx9DhgzhwQcfpEoVbYhU9u9Xa+LGGAfy5zsfS/6CJVeBqkAqsAX4QESOlXKcNloTV0oV18mTJ3n22WcJDw8nJycHDw8PBg4cyKxZs2wLuihVXt3rLWY7gGbAi0ADEfEWkXpAMLAbmGuMGVdi0SqlVAlr2rQp//73v0lLS2PdunUMGjSIzZs32yaV+e6771i8eDGJiYkWR6pU8RSlJu4kIoXemFmUY0qK1sSVUiXh2rVrODo6Yoxh5syZzJ8/H4BOnToxfPhwRowYQcuWha6MrFSZuKeaeFGSc1klcKWUKilVqlTBGAPAW2+9xeHDh5k3bx5OTk68/PLL9O3bl18qOUlJSdjbIGBVORR5dLoxJgh4GfAhf0CcAURE2pZeeLfSmrhSqrSdOnWKuLg4unfvTm5uLg0bNqR69eq2W9weeOABHBx0wktVNkpkshdjzBFgBhAN2FYnEJGEkgiyqDSJK6XKUnZ2NqtXr+bzzz9n69atZGdn06BBA95++20ef/xxq8NTlUBJzZ2eKiIbRSRORBJ+eZRQjEopVS45Ozvz1FNP8eWXX5KSksInn3xC9+7dbdO97tu3j/Hjx7N582ays7MtjlZVNsWpiYcAjwPbyb/NDAAR+bx0Qrs9rYkrpcqTzz77jClTpnDhwgXc3d0ZNmwYo0aNIjQ0VO9FVyWipJrTVwOtgEP8X3O6iMhTJRJlEWkSV0qVN9nZ2WzdupWwsDC++OILRISUlBRcXFyIj4+ncePGmtDVXSupJB4tIoElGtld0CSulCrPsrKy+PHHH+nYsSMA/v7+pKamMmLECEaNGkWPHj10kRZVLCXVJ77bGKNTGymlVCFcXFxsCVxE+Mtf/kJISAgrV66kd+/eeHl5sWzZMoujVBVFcZJ4dyDKGHPEGHPQGBNtjDlYWoEppZS9M8YwbNgwPv30U1JTU1m7di3BwcG2ldcSExOZPn06kZGReh+6uivFaU73ud12vcVMKaXuzvr163n88cfJycmhZcuWjBkzhjFjxnDfffdZHZoqR+6pT9wYY+RXDirKMSVFk7hSqiL5+eefWbduHWvWrOGbb77B0dGRlJQUPDw8yM3N1f5zde8LoBhjphpjmtxUqLMxprcx5mNgfEkEqpRSlU3t2rWZPHkyERERJCYm8sknn+Dh4QHAgAED6NevHytXruTSpUsWR6rKo6Ik8f5ALvCpMea0MeZHY8xJ4Bj5942/IyIrSjFGpZSqFLy9vRk1ahSQPyiuS5cuHD16lPHjx1O/fn0ee+wxvv32W4ujVOVJURZAyRKRJSLSjfx500OAjiLiIyKTRSSqtINUSqnKxhjDnDlzOHnyJLt27WLixIls27aNvXv3AnD58mX27t2rA+IquSIPbCsvtE9cKVVZ5eTkkJOTg6urKx9//DETJkygdevWjB8/nnHjxuHl5WV1iKoUlNR94koppSzk5OSEq6srAEOHDuWDDz7Aw8ODWbNm0aRJE/r168fly5ctjlKVJU3iSillh9zc3Jg8eTK7du3i2LFjvPzyy1SrVo0aNWoA8NFHH/Htt99qc3sFV5z7xHsDY4ELQAxwEIgRkauFva6kaXO6UkoVLicnB29vb5KTk/Hz8+PJJ5/kySefpGnTplaHpu5CSTWnrwa+BHYDTYE/k78YSlEC8DbG7DDGxBpjDhlj/lCw3cMYs9UYc6zg39rFiEcppdRtODk5ceLECVauXEnTpk2ZM2cOzZo147333rM6NFXCilMT3ykiPe7qJMY0BBqKyH5jTE1gHzAUmACcF5G5xphZQG0ReaGwsrQmrpRSxZOYmMjq1asZNmwYrVu3ZufOnaxZs4ZJkyYRFBSEMcbqEFUh7qkmboxZaYx5HvivMeaPdxOAiJwRkf0FP18CYgEv4BHg44LDPiY/sSullCpBTZo04aWXXqJ169YAxMbGsmrVKu6//37atWvHwoULSUtLszhKdTeK0pz+ccFxDYAnjDEJxpiNxpjXjTEji3tCY4wv0AHYA9QXkTOQn+iBesUtTymlVPE888wznDlzhqVLl1K1alX+8Ic/EBQUpIPg7FCx7xM3xlQB2gDtgLYiMqMYr60BfAO8ISKfG2MuiIj7dft/FpFb+sWNMU8DTwM0adKkU0JCma65opRSFdr//u//kpiYyODBg8nNzaVXr16EhoYyceJEvL29rQ6v0runBVBKMAgn8gfG/UdE3i7YdgToKSJnCvrNI0SkZWHlaJ+4UkqVnuTkZMaNG8e2bdswxtCvXz8mTZrE4MGDcXZ2tjq8SsnyyV5M/qiJj4DYXxJ4gY383+Ip44F/lUU8Simlbq9+/fps3bqVkydP8qc//Yno6GgeffRRwsPDAbTJvZwpk5q4MaY78C0QDeQVbH6J/H7xMKAJkAiMFJHzhZWlNXGllCo7ubm5bNu2jdDQUBwdHXnppZc4cOAAU6ZMYeDAgVSpUsXqECu8wmriZXL1RWQXcKd7GELKIgallFLF5+joSL9+/WzP69Wrx8GDBxk6dCheXl5MnjyZSZMm6bztFtFpV5VSShXZ888/T0JCAhs2bCAgIIDZs2fzP//zP7b9eXl5hbxalTRN4koppYqlSpUqDB06lH//+9+cOHGCOXPmABATE0OzZs148803SU5OtjjKykGTuFJKqbvWtGlTmjdvDsDVq1fx8/PjpZdeonHjxowaNYqIiAgdDFeKNIkrpZQqEZ06dSI8PJzY2FimTp3Ktm3bGDhwIBcvXgR0ZHtp0CSulFKqRLVq1Yq3336bU6dO8fXXX+Pm5oaIEBISwrRp0zh27JjVIVYYmsSVUkqVimrVqtGtWzcAMjMzadiwIYsXL6ZFixb079+fL7/8ktzcXIujtG+axJVSSpU6V1dX1qxZQ2JiInPmzCE6OprBgwfz2WefWR2aXdMkrpRSqsw0aNCAV155hfj4eMLCwhg2bBgAS5cuZdKkSURFRVkboJ3RJK6UUqrMOTk5MXLkSKpVqwZAamoqn3zyCR06dKB79+6EhYVx7do1i6Ms/zSJK6WUstwrr7zCqVOnePvttzl79iyjR49m7NixVodV7mkSV0opVS7Url2badOmceTIETZu3MjUqVMBOHPmDFOnTuXo0aMWR1j+aBJXSilVrjg6OjJ48GC6d+8OwH//+18++OADWrZsyeDBg9m+fbvec15Ak7hSSqlybfjw4SQmJvLqq6+yZ88eQkND6dixI1evXrU6NMuVyVKkJel2S5Hm5OSQlJREVlaWRVEpVf64uLjQuHFjnJycrA5FqRKTlZXFp59+ypEjR5g7dy4Aq1atok+fPjRo0MDi6EpHYUuRVogkHhcXR82aNalTpw7G3GnFU6UqDxEhLS2NS5cu4efnZ3U4SpWa06dP4+3tjaOjI48//jjTpk2jffv2VodVogpL4hWiOT0rK0sTuFLXMcZQp04dbZ1SFV6jRo04fPgwU6ZMYf369XTo0IHQ0NBKMwiuQiRxQBO4UjfRz4SqLJo3b87ChQtJSkpi3rx5JCQk4OHhAUBiYiLZ2dkWR1h6KkwSt9rChQtp3bp1ubivcfbs2SxYsOCeyoiIiGDQoEGFHnPhwgWWLFlie3769GkeffTRezpvWfnzn//Mtm3bSvUcRfl/ePfdd7ly5Yrt+YABA7hw4UKxz7VixQqee+65Yr9OqYrE3d2dmTNncvToUTw9PRERRo4cSdOmTVmwYIFtNbWKRJN4CVmyZAlbtmxhzZo1RTq+pGYiEhHy8vJKpKziujmJN2rUiHXr1lkSS3HNmTOH0NBQq8O4JYlv2bIFd3d36wJSqgK4vhVq9uzZtGjRghkzZuDt7c3MmTNJSkqyMLqSpUm8BEyZMoWTJ08yZMgQ3nnnHc6fP8/QoUNp27YtXbt25eDBg0D+L9PTTz9N3759efLJJxkwYIBtX4cOHZgzZw6QP3PRhx9+yOXLlwkJCaFjx44EBgbyr3/9C4D4+Hhat27Ns88+S8eOHfnpp5944403aNmyJaGhoRw5cuS2ca5du5aAgADatWtHjx49gPzxBBMnTiQwMJAOHTqwY8eOW153c40yICCA+Ph4Zs2axYkTJ2jfvj0zZswgPj6egICAQstdsWIFw4cPp3///jRv3pyZM2f+6vVduHAhbdq0oW3btjz22GMAZGRk8NRTT9G5c2c6dOhguzYrVqxg6NChDB48GD8/P9577z3efvttOnToQNeuXTl//jwAEyZMuO0Xjl+75pMnT8bf35++ffuSmZkJwLJly+jcuTPt2rVjxIgRNyRlgBMnTtCxY0fb82PHjtGpUycWLlzI6dOn6dWrF7169QLA19eXc+fOAbBy5Uratm1Lu3bteOKJJwDYtGkTXbp0sfX7JScn/+r1U6qyMsbw8MMPEx4ezt69exkwYAB/+9vf2LRpE1BB1jcXEbt6dOrUSW72448/3vD8oYceuuWxePFiERHJyMi47f7ly5eLiEhqauot+4rCx8dHUlNTRUTkueeek9mzZ4uIyPbt26Vdu3YiIvLqq69Kx44d5cqVKyIi8uabb8p7770n6enpEhQUJH379hURkZ49e8rhw4clJydH0tPTbXE1a9ZM8vLyJC4uTowx8t///ldERCIjIyUgIEAyMjIkPT1dmjVrJvPnz78lxoCAAElKShIRkZ9//llERBYsWCATJkwQEZHY2Fjx9vaWzMxM2bFjhwwcONAW9/Xl+fv7S1xcnMTFxYm/v79t+/XP71Tu8uXLxc/PTy5cuCCZmZnSpEkTSUxMLPTaNmzYULKysm6I+8UXX5RVq1bZtjVv3lwuX74sy5cvl2bNmsnFixclJSVFatWqJe+//76IiDz//PPyzjvviIjI+PHjZe3atbecq7Br7ujoKAcOHBARkZEjR9rOf+7cOdvrX375ZVm4cOEt161nz56217744ou2Y67/vbn+eUxMjLRo0cK2Ly0tTUREzp8/L3l5eSIismzZMpk+fbqIiCxfvlx+97vf3fJ+bv5sKFXZxcXFSUZGhoiILFmyRPr37y/bt2+3fa7KIyBS7pATtSZeCnbt2mWrOfXu3Zu0tDTS09MBGDJkiG3C/+DgYHbu3MmuXbsYOHAgly9f5sqVK8THx9OyZUtEhJdeeom2bdsSGhrKqVOnbDUvHx8funbtCsC3337LsGHDcHV1pVatWgwZMuS2cXXr1o0JEyawbNky2xq+18faqlUrfHx8SmRUZ2HlhoSE4ObmhouLC23atCEhIaHQstq2bcvYsWNZvXo1VapUAeDrr79m7ty5tG/fnp49e5KVlUViYiIAvXr1ombNmtStWxc3NzcGDx4MQGBgIPHx8YWeq7Br7ufnZ7t1pVOnTrayYmJiCA4OJjAwkDVr1nDo0KFbyp00aRLLly8nNzeXzz77jDFjxhQaR3h4OI8++iienp4AtkE6SUlJ9OvXj8DAQObPn3/bcyml7szX1xdXV1cgfxGW/fv3ExISQlBQEGvXrrW79c2rWB1AaYiIiLjjPldX10L3e3p6Frq/KOQ2TTS/9NFUr17dtq1z585ERkbStGlT+vTpw7lz51i2bBmdOnUCYM2aNaSmprJv3z6cnJzw9fW13TJ0fTnXl1+YpUuXsmfPHjZv3kz79u2JiooqUnNSlSpVbuh3L8ptS4WVW7VqVdvPjo6Ovzo+YPPmzezcuZONGzfy+uuvc+jQIUSE9evX07JlyxuO3bNnzw3lOzg42J47ODjccq49e/bwzDPPAPn95OfPn7/jNb857l+a0ydMmMAXX3xBu3btWLFixW1/f0aMGMFrr71G79696dSpE3Xq1Cn0PYvIbf9Pp06dyvTp0xkyZAgRERHMnj270HKUUnc2adIkxo0bx6pVq5g/fz6jRo1i9OjR/POf/7Q6tCLTmngp6NGjh22AW0REBJ6entSqVeuW45ydnfH29iYsLIyuXbsSHBzMggULCA4OBiA9PZ169erh5OTEjh077lhj7dGjBxs2bCAzM5NLly7Z+ntuduLECbp06cKcOXPw9PTkp59+uiHWo0ePkpiYeEti9PX1Zf/+/QDs37+fuLg4AGrWrMmlS5d+9RrcqdybPfnkk/zwww83bMvLy+Onn36iV69evPXWW1y4cIHLly/Tr18/Fi1aZPuycODAgULLvpMuXboQFRVFVFQUQ4YMKfI1v96lS5do2LAhOTk5dxzY6OLiQr9+/fjtb3/LxIkTbdvvdA1DQkIICwsjLS0NwNaXn56ejpeXFwAff/xxsd+vUupGLi4uTJ48mdjYWMLCwmxf6s+dO8e7777L5cuXLY6wcJrES8Hs2bOJjIykbdu2zJo1q9A/tsHBwdSvXx9XV1eCg4NJSkqyJfGxY8cSGRlJUFAQa9asoVWrVrcto2PHjowePZr27dszYsQI2+tvNmPGDAIDAwkICKBHjx60a9eOZ599ltzcXAIDAxk9ejQrVqy4ocYJ+bXI8+fP0759e95//31atGgBQJ06dejWrRsBAQHMmDHjhtcUpdybHTx4kIYNG96wLTc3l3HjxtkGyE2bNg13d3deeeUVcnJyaNu2LQEBAbzyyiuFll1URb3m13v99dfp0qULffr0KfT4sWPHYoyhb9++tm1PP/00Dz/8sG1g2y/8/f15+eWXeeihh2jXrh3Tp08H8n+3Ro4cSXBwsK2pXSl17xwdHRk5cqTts7hhwwamTZuGj48Pr732mu2LdHlTIaZdjY2NpXXr1hZFpErCxYsX+c1vfsPatWutDqXULFiwgPT0dF5//fUyO6d+NpS6e99//z1z585l06ZNVK9end/+9rfMmzcPB4eyrf8WNu1qhewTV/anVq1aFTqBDxs2jBMnThAeHm51KEqpInrwwQfZuHEj0dHRzJs3j7i4OFsCT0lJoV69ehZHqElcqTKxYcMGq0NQSt2lwMBAVq9ebRvge+zYMfz9/RkxYgSzZs2iXbt2lsWmfeJKKaVUEfxSC3d3d2fatGm2O30GDRrE7t27rYnJkrOWAnvr21eqtOlnQqnSUbduXebNm0diYiKvv/46u3fvpmfPnra7ScpShUjiLi4upKWl6R8tpQpIwXriLi4uVoeiVIXl7u7On/70J+Lj4/nyyy9/df6H0mB5n7gxpj/wd8AR+FBE5ha3jMaNG5OUlERqamqJx6eUvXJxcaFx48ZWh6FUhVejRg3LFlSyNIkbYxyBxUAfIAnYa4zZKCI/FqccJycn/Pz8SiNEpZRSqtyyujn9fuC4iJwUkWzgn8AjFseklFJK2QWrk7gX8NN1z5MKtimllFLqV1idxG+3ascto9OMMU8bYyKNMZHa762UUkrls3pgWxLgfd3zxsDpmw8SkQ+ADwCMManGmF9flaLoPIFzJVheZaTXsGTodbx3eg3vnV7De1fS19DnTjssnTvdGFMFOAqEAKeAvcAYESmzRZKNMZF3mpNWFY1ew5Kh1/He6TW8d3oN711ZXkNLa+Iics0Y8xzwH/JvMft/ZZnAlVJKKXtmdXM6IrIF2GJ1HEoppZS9sXpgW3nwgdUBVAB6DUuGXsd7p9fw3uk1vHdldg3tbj1xpZRSSuXTmrhSSillpyp1EjfG9DfGHDHGHDfGzLI6HntjjPE2xuwwxsQaYw4ZY/5gdUz2yhjjaIw5YIz50upY7JExxt0Ys84Yc7jg9/EBq2OyN8aYaQWf4xhjzKfGGF09pwiMMf/PGJNijIm5bpuHMWarMeZYwb+1S+v8lTaJXzdv+8NAG+BxY0wba6OyO9eAP4pIa6Ar8Du9hnftD0Cs1UHYsb8D/xaRVkA79FoWizHGC/g9ECQiAeTfLfSYtVHZjRVA/5u2zQK2i0hzYHvB81JRaZM4Om/7PRORMyKyv+DnS+T/4dRpc4vJGNMYGAh8aHUs9sgYUwvoAXwEICLZInLB0qDsUxWgWsH8Ha7cZuItdSsR2Qmcv2nzI8DHBT9/DAwtrfNX5iSu87aXIGOML9AB2GNxKPboXWAmkGdxHPaqKZAKLC/okvjQGFPd6qDsiYicAhYAicAZIF1EvrY2KrtWX0TOQH5lB6hXWieqzEm8SPO2q19njKkBrAeeF5GLVsdjT4wxg4AUEdlndSx2rArQEXhfRDoAGZRi82VFVNBn+wjgBzQCqhtjxlkblSqKypzEizRvuyqcMcaJ/AS+RkQ+tzoeO9QNGGKMiSe/S6e3MWa1tSHZnSQgSUR+aQVaR35SV0UXCsSJSKqI5ACfAw9aHJM9SzbGNAQo+DeltE5UmZP4XqC5McbPGONM/iCOjRbHZFeMMYb8fshYEXnb6njskYi8KCKNRcSX/N/BcBHRGlAxiMhZ4CdjTMuCTSHAjxaGZI8Sga7GGNeCz3UIOjjwXmwExhf8PB74V2mdyPJpV62i87aXiG7AE0C0MSaqYNtLBVPpKlWWpgJrCr6QnwQmWhyPXRGRPcaYdcB+8u86OYDO3FYkxphPgZ6ApzEmCXgVmAuEGWN+Q/4XpJGldn6dsU0ppZSyT5W5OV0ppZSya5rElVJKKTulSVwppZSyU5rElVJKKTulSVwppZSyU5rElVJKKTulSVwppZSyU5rElVI3MMZUM8Z8U7DGue/16yQXsxxnY8zOglWxlFKlQJO4UupmTwGfi0juvRRSsMTvdmB0iUSllLqFJnGlKhFjzA5jTJ+Cn/9ijFl4m8PGcpu5no0xTQuW+uxcUEM/XLDsZ4wxZo0xJtQY850x5pgx5v6Cl31RUJ5SqhRoEleqcnkVeNkYM5b89d+nXb+zYO7xpiISf9P2luSvVjdRRPYWbL4P+DvQFmgFjAG6A/8DvFRwTAzQuVTeiVKq8i6AolRlJCI7C1apmg70vE2TuSdw4aZtdcmvmY+4aZGgOBGJBjDGHAK2i4gYY6IB34Lz5Rpjso0xNUXkUsm/I6UqN62JK1WJGGMCgYbA1Tsk1UzA5aZt6cBP5K9ad72r1/2cd93zPG6sIFQFsu42ZqXUnWkSV6qSMMY0BNYAjwAZxph+Nx8jIj8DjsaY6xN5NjAUeNIYM6aY56wDpIpIzl0HrpS6I03iSlUCxhhX4HPgjyISC7wOzL7D4V+T37dtIyIZwCBgmjHmkWKcuheg68srVUp0PXGl1A2MMR2A6SLyRAmU9TnwoogcuffIlFI305q4UuoGInIA2GGMcbyXcgpGun+hCVyp0qM1caWUUspOaU1cKaWUslOaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslP/H82/ICYOVs+GAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from os import getcwd, chdir\n",
    "cwd = getcwd()\n",
    "chdir(r'../../timml')\n",
    "import timml\n",
    "chdir(cwd)\n",
    "\n",
    "dpar = 1e-4\n",
    "M0 = timml.Model3D(kaq=K, z=[0., -b])\n",
    "timml.ImpLineDoublet1D(M0, xld=0.)\n",
    "timml.HeadLineSink1D(M0, xls=L, hls=BC3h, res=BC3c/100., wh=1.)\n",
    "timml.StripAreaSink(M0, 0., L, R)\n",
    "M0.solve(silent=True)\n",
    "H0 = M0.headalongline(X, 0.).flatten()\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., H0, 'k--', mfc='none', label='forward solution, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel('$h$ (m)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAChCAYAAABakLybAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAh30lEQVR4nO3deXxV5Z3H8c+PRZBFRUBAsIRNQxKSAJFFlOCCIFIRBZWhSLBIddDqWGkpdgbnZfsaOzBuYwtSi6iDC1pRWmm1RRBo3cIiAoKABAggm7JEpBL4zR/35HrJzUqSG8L9vl+v+7r3nN9zznnO4YT7u88553nM3REREZH4Vau6KyAiIiLVS8mAiIhInFMyICIiEueUDIiIiMQ5JQMiIiJxTsmAiIhInItpMmBmM81st5mtLiZuZvaEmW00s1Vm1i0ilmNmn5jZSjPLjph/rpn91cw2BO9NYrEvIiIip4tYtwzMAgaWEL8G6BS8xgHTCsUvd/d0d8+ImDcRWODunYAFwbSIiIiUUUyTAXdfDHxZQpEhwHMe8j5wjpm1KmW1Q4Bng8/PAtdXuKIiIiJx5FS7Z6A1sC1iOjeYB+DA22a2zMzGRZRp4e47AYL382JSUxERkdNEnequQCFWxLyC/pL7uPsOMzsP+KuZrQtaGsq+8lASMQ6gYcOG3RMTEytWWxERkRpi2bJle929eVGxUy0ZyAUuiJhuA+wAcPeC991mNhfoASwGdplZK3ffGVxS2F3cyt19BjADICMjw7Ozs4srKiIicloxsy3FxU61ywTzgFuDpwp6AQeCL/mGZtYYwMwaAlcDqyOWGR18Hg28EetKi4iI1GQxbRkwsxeBfkAzM8sFJgN1Adx9OjAfGARsBA4DY4JFWwBzzaygzi+4+1+C2MPAHDP7IbAVGB6TnRERETlNxDQZcPcRpcQdGF/E/M+BtGKW2QdcWSkVFBERiUOn2j0DInIKO3r0KLm5uRw5cqS6qyIixahfvz5t2rShbt26ZV5GyYCIlFlubi6NGzcmISGB4LKdiJxC3J19+/aRm5tLu3btyrzcqXYDoYicwo4cOULTpk2VCIicosyMpk2blrv1TsmAiJSLEgGRU9vJ/I0qGRCRuPL666+zdu3aci+3aNEiBg8eXK5lHnvsMQ4fPhyeHjRoEPv37y/3tmNp+vTpPPfccwDMmjWLHTt2hGNjx44t9dhdcsklAOTk5PDCCy9UXUVPwo4dOxg2bFiVb6dRo0YlxleuXMn8+fPD0/PmzePhhx8+qW0lJCSwd+/ek1o2kpIBEYkb+fn5J5UM5Ofnn9T2CicD8+fP55xzzjmpdcXKHXfcwa233gpEJwNPP/00SUlJJS7/j3/8Azg1k4Hzzz+fV199tbqrEZUMXHfddUycWL1j7CkZEJEaIycnh8TEREaPHk1qairDhg0Lf9kuW7aMzMxMunfvzoABA9i5cycA/fr1Y9KkSWRmZvLrX/+aefPmMWHCBNLT09m0aRP9+vWjoDfSvXv3kpCQAIS+CIcPH873v/99rr76agAOHjzI0KFDSUpK4o477uD48eMA3HnnnWRkZJCcnMzkyZMBeOKJJ9ixYweXX345l19+OXDir7hHHnmElJQUUlJSeOyxx8L717lzZ26//XaSk5O5+uqr+eabb6KOwyuvvEJKSgppaWn07dsXgGPHjjFhwgQuvvhiUlNTeeqpp4BQi0a/fv0YNmwYiYmJjBw5ktBT3DBx4kSSkpJITU3l/vvvB+DBBx9k6tSpvPrqq2RnZzNy5EjS09P55ptvwsdq2rRp/PSnPw3XZ9asWdx9993Ad7+KJ06cyJIlS0hPT+fRRx/lsssuY+XKleFl+vTpw6pVq4r9t965cyd9+/YlPT2dlJQUlixZAsDbb79N79696datG8OHDycvLy98bCdNmkTv3r3JyMhg+fLlDBgwgA4dOjB9+vTw8U1JSSlye7/73e+4+OKLSUtL48YbbwyfV1lZWfz4xz/mkksuoX379uFkIi8vjyuvvJJu3brRpUsX3ngjur+7UaNGnTB/5MiRzJs3j//4j//g5ZdfJj09nZdffplZs2Zx1113AbBr1y6GDh1KWloaaWlp4eTq+uuvp3v37iQnJzNjxoxij9tJc/e4fHXv3t1FpHzWrl1brdvfvHmzA7506VJ3dx8zZoxPmTLFv/32W+/du7fv3r3b3d1feuklHzNmjLu7Z2Zm+p133hlex+jRo/2VV14JT2dmZvpHH33k7u579uzxtm3burv7M888461bt/Z9+/a5u/vChQu9Xr16vmnTJs/Pz/errroqvJ6CMvn5+Z6Zmekff/yxu7u3bdvW9+zZE95WwXR2dranpKR4Xl6eHzp0yJOSknz58uW+efNmr127tq9YscLd3YcPH+7PP/981HFISUnx3Nxcd3f/6quv3N39qaee8oceesjd3Y8cOeLdu3f3zz//3BcuXOhnnXWWb9u2zY8dO+a9evXyJUuW+L59+/zCCy/048ePn7CeyZMn+5QpU6KOTeT07t27vUOHDuH5AwcO9CVLlri7e8OGDcPH69prrw2XmTVrlt9zzz3u7r5+/Xov7f/gqVOn+i9/+cvwcT148KDv2bPHL7vsMs/Ly3N394cfftj/8z//M3xsf/vb37q7+7333utdunTxgwcP+u7du7158+buHjp/kpOTi9ze3r17w58feOABf+KJJ9w9dL4MGzbMjx075mvWrAnv99GjR/3AgQPuHjpvOnToED6WBcdg0aJFPmTIEHd3379/vyckJPjRo0f9mWee8fHjx4e3Fzl90003+aOPPhre7/3797v7d+fY4cOHPTk5OVzfwudYgaL+VoFsL+Y7UY8WishJ69evX9S8m266iX/913/l8OHDDBo0KCqelZVFVlYWe/fujbp+u2jRolK3ecEFF9CnTx8AfvCDH/DEE08wcOBAVq9eTf/+/YHQr+RWrb4b/fzmm28ux159p3///px77rnh6R49etC+fXsARowYwdKlSxk2bBhz5sxhxowZ5Ofns3PnTtauXUtqamqx6126dClDhw6lYcOGANxwww0sWbKE6667jnbt2pGeng5A9+7dycnJiVq+T58+ZGVlcdNNN3HDDTcAoV/Mq1atCv9yPXDgABs2bOCMM86gR48etGnTBoD09HRycnLo1asX9evXZ+zYsVx77bXluh+iefPmtG/fnvfff59OnTqxfv368L9JcYYPH85DDz3ElClTmDlzJllZWSWWv/jii7nttts4evQo119/Penp6bz77rusXbs2vK1vv/2W3r17h5e57rrrAOjSpQt5eXk0btyYxo0bU79+/VLv1Vi9ejW/+MUv2L9/P3l5eQwYMCAcu/7666lVqxZJSUns2rULCP2QnjRpEosXL6ZWrVps376dXbt20bJly/BymZmZjB8/nt27d/Paa69x4403UqdOyV+777zzTviejdq1a3P22WcDoZamuXPnArBt2zY2bNhA06ZNS1xXeSgZEJEapfCd0maGu5OcnMx7771X5DIFX7pFqVOnTri5v/DjWIWXK2rbmzdvZurUqXz00Uc0adKErKysUh/rcvdiY/Xq1Qt/rl27dpGXCaZPn84HH3zAm2++SXp6OitXrsTd+d///d8TvsQglGAVXmd+fj516tThww8/ZMGCBbz00ks8+eSTvPPOOyXWO9LNN9/MnDlzSExMZOjQoaXewd6gQQP69+/PG2+8wZw5cyhtoLi+ffuyePFi3nzzTUaNGsWECRNo0qQJ/fv358UXXyxymYL9rFWr1gn7XKtWraj7PsaMGcOKFSs4//zzmT9/PllZWbz++uukpaUxa9asExLTyHUV/NvNnj2bPXv2sGzZMurWrUtCQkKR/+6jRo1i9uzZvPTSS8ycObPEfS7OokWL+Nvf/sZ7771HgwYN6NevX6V3/KVkQEROWkm/5Bs0aFBivFmzZmVqCShs69atvPfee/Tu3ZsXX3yRSy+9lIsuuog9e/aE5x89epTPPvuM5OTkqOUbN27MoUOHwtMJCQksW7aMHj16lHpz2YcffsjmzZtp27YtL7/8MuPGjePgwYM0bNiQs88+m127dvHnP/853GJSsK1mzZqdsJ6+ffuSlZXFxIkTcXfmzp3L888/X+ZjsGnTJnr27EnPnj354x//yLZt2xgwYADTpk3jiiuuoG7dunz22We0bt262HXk5eWFW2969epFx44dSz1WkW644QZ+9atf0bZtW37961+XadmxY8fy/e9/n8suuyzc4vLhhx/y5JNPhn8NF9iyZQutW7fm9ttv5+uvv2b58uU88MADjB8/no0bN9KxY0cOHz5Mbm4uF154YanHrLBnnnnmhOlDhw7RqlUrjh49yuzZs0s8dhBqeTnvvPOoW7cuCxcuZMuWogcEzMrKokePHrRs2TJ8PpZ0XK+88kqmTZvGvffey7Fjx/j66685cOAATZo0oUGDBqxbt47333+/3PtbGt1AKCI1SufOnXn22WdJTU3lyy+/5M477+SMM87g1Vdf5Wc/+xlpaWmkp6eHb7wq7JZbbmHKlCl07dqVTZs2cf/99zNt2jQuueSSUh/R6t27NxMnTiQlJYV27dqFb/Tq2rUrycnJ3HbbbSc0l48bN45rrrkmfANhgW7duoW/JHr27MnYsWPp2rVrmY/BhAkT6NKlCykpKfTt25e0tDTGjh1LUlIS3bp1IyUlhR/96EclPgVx6NAhBg8eTGpqKpmZmTz66KNRZbKysrjjjjvCNxBGatKkCUlJSWzZsoUePXpELZuamkqdOnVIS0sLr7t79+6cddZZjBkzJlxu69atnHnmmVHLL1q0iPT0dLp27cof/vAH7rnnHpo3b86sWbMYMWIEqamp9OrVi3Xr1pX5uJXkoYceomfPnvTv35/ExMRSy48cOZLs7GwyMjKYPXt2scu0aNGCzp07n7DPl19+OWvXrg3fQBjp8ccfZ+HChXTp0oXu3buzZs0aBg4cSH5+Pqmpqfz7v/87vXr1qtjOFsFKaq46nWVkZHhpzVQicqJPP/2Uzp07V9v2c3JyGDx4MKtXry69sJxyduzYQb9+/Vi3bh21aoV+i06YMIFRo0aVeI9FTXb48GG6dOnC8uXLw9f/Y6Gov1UzW+buGUWVV8uAiIhUueeee46ePXvyq1/9KpwIAEyZMuW0TQT+9re/kZiYyN133x3TROBkqGVARMqsulsGRKRs1DIgIiIi5aJkQETKJV5bE0VqipP5G1UyICJlVr9+ffbt26eEQOQU5e7s27eP+vXrl2u5mPYzYGYzgcHAbneP6iDaQr1WPA4MAg4DWe6+3MwuAJ4DWgLHgRnu/niwzIPA7cCeYDWT3H1+4XWLSMW1adOG3Nxc9uzZU3phEakW9evXD/c4WVax7nRoFvAkoS/2olwDdApePYFpwXs+8JMgMWgMLDOzv7p7wdBjj7r71CqtuYhQt25d2rVrV93VEJFKFtPLBO6+GPiyhCJDgOeCMRXeB84xs1buvtPdlwfrOAR8CpTcPZSIiIiUyal2z0BrYFvEdC6FvvTNLAHoCnwQMfsuM1tlZjPNrEmV11JEROQ0cqolA0WNdBG+U8nMGgF/AO5194PB7GlAByAd2An8T7ErNxtnZtlmlq1rniIiIiGnWjKQC1wQMd0G2AFgZnUJJQKz3f21ggLuvsvdj7n7ceB3QHQn2d+VneHuGe6e0bx58yrZARERkZrmVEsG5gG3Wkgv4IC77wyeMvg98Km7PxK5gJm1ipgcCqjTchERkXKI9aOFLwL9gGZmlgtMBuoCuPt0YD6hxwo3Enq0sGCYpz7AKOATM1sZzCt4hPC/zSyd0OWEHOBHMdgVERGR00ZMkwF3H1FK3IHxRcxfStH3E+DuoyqndiIiIvHpVLtMICIiIjGmZEBERCTOKRkQERGJc0oGRERE4pySARERkTinZEBERCTOKRkQERGJc5WSDJjZ1ZWxHhEREYm9k+50KOhNcAWhzoAuBd6urEqJiIhI7FSkB8LfuvsSADP7cyXVR0RERGKsXJcJzOz+gs8FiUDweVVlVkpERERip7z3DLSvklqIiIhItSnvZYK2ZjYockYwcqCIiIjUUOVNBuYAzfhuBEGv3OqIiIhIrJUpGTCzK4CRwH5gNbAKWO3u/6y6qomIiEgslLVl4P+A8UH5VOB6IBnoWDXVEhERkVgpazKw0d3nBp9fqarKiIiISOyV+DSBmT1nZvcC75nZT2JTJREREYml0h4tfDYo0xIYZWZbzGyemT1kZsOrvnoiIiJS1UpMBtx9gbs/4u6j3T0d6AD8AvgM6FHejZnZTDPbbWari4mbmT1hZhvNbJWZdYuIDTSz9UFsYsT8c83sr2a2IXhvUt56iYiIxLNydTrk7vnuvsrdn3f3CSexvVnAwBLi1wCdgtc4YBqAmdUGfhPEk4ARZpYULDMRWODunYAFwbSIiIiUUak3EJpZPeCnQG/ga2A5MMfdN5V3Y+6+2MwSSigyBHjO3R1438zOMbNWQAKhmxg/D+r0UlB2bfDeL1j+WWAR8LPy1q2i3n//fXbt2nXCvIYNG3LVVVcBsHTpUvbt23dC/Oyzz6Zfv34ALFq0iAMHDpwQb9q0KZdeeikACxYsIC8v74R4ixYt6NWrFwBvvfUWR44cOSHeunVrMjIyAHjzzTfJz88/Id62bVvS09MBeOONN6L2qUOHDqSkpHD06FHmz4/uW+qiiy4iMTGRI0eO8NZbb0XFk5OT6dixI3l5eSxYsCAqnpaWRkJCAvv37+fdd9+Ninfv3p02bdqwd+9e/v73v0fFe/bsScuWLfniiy/44IMPouJ9+vShWbNm5ObmsmzZsqh4ZmYm55xzDjk5OXz88cdR8SuuuILGjRuzceNG1qxZExW/+uqrOfPMM1m3bh3r16+Pig8aNIi6deuyevVqNm2K/nMZMmQIACtXrmTLli0nxOrUqcO1114LQHZ2Ntu3bz8hXq9ePQYODOXVOvd07hWmc+/0OPcaNGhAgwYNospUCXcv8QU8ATwNpAGPAlOBxcAPSlu2mPUlEOqjoKjYn4BLI6YXABnAMODpiPmjgCeDz/sLreOrErY9DsgGsr/3ve95ZRo8eLAT6oQp/OrUqVM4npmZGRXv1q1bON6tW7eoeGZmZjh+4YUXRsUHDx4cjrdq1Soqfsstt4TjjRo1iorffvvt4XjhGOD33Xefu7sfPHiwyPjkyZPd3X3Hjh1FxqdOneru7uvXry8yPn36dHd3z87OLjI+e/Zsd3dfuHBhkfF58+a5u/u8efOKjC9cuNDd3V944YUi49nZ2e7u/tRTTxUZX7dunbu7T506tcj49u3b3d198uTJRcYPHjzo7u733XdfkfEC48aNi4o1atQoHB8xYkRUvGXLljr3dO7p3DvNz70tW7Z4ZQKyvZjvRwvFi2dm7wJD3H2/mf3d3fsErQWL3L13iQsXvb4E4E/unlJE7E3gv9x9aTC9gFCrRHtggLuPDeaPAnq4+91mtt/dz4lYx1fuXup9AxkZGZ6dnV3e6hdr06ZNHDx48IR59erVIykpdDVjw4YNURnumWeeSWJiIgDr1q3jm2++OSHeqFEjOnXqBMDatWv55z9P7OPprLPOokOHDgB88sknURlwkyZNSEhIAODjjz/m+PHjJ8SbNm3K9773PQBWrFgRtU/nnXcerVu35tixY6xaFT0WVcuWLWnVqhVHjx5l9ero20DOP/98WrRowZEjR/j000+j4hdccAHNmjXj8OHDRf66adu2Leeeey6HDh1i48aNUfH27dtz9tlnc+DAAT7//POoeMeOHWncuDFffvll1K8fCGX4DRo0YO/evWzbti0q3rlzZ+rXr8+uXbvYsWNHVDwlJYW6deuyc+dOvvjii6h4amoqtWvXZvv27ezevTsq3rVrVwC2bt0a9eupVq1apKWlAZCTk8NXX311QrxOnTp06dIF0Lmnc0/nXqTT6dyrV68eZ5xxRlSZk2Vmy9w9o8hYGZKBwcA/3P1LM3sD2ASsA0a7e5+TqEwCxScDTxFKMl4MptcTugSQADzo7gOC+T8HcPf/Kijj7juDSwqL3P2i0upR2cmAiIjIqaykZKC0fgZ6A2+6+5fBrOGEmtkbBZ8r2zzg1uCpgl7AAXffCXwEdDKzdmZ2BnBLULZgmdHB59FA9EUgERERKVZpNxCOBn5jZp8BfwH+4u4vnOzGzOxFQr/0m5lZLjAZqAvg7tOB+cAgYCNwGBgTxPLN7C7gLaA2MNPdC+6qeRiYY2Y/BLZSNUmKiIjIaavUywQAZpZI6LG+AcDZwEJCycHf3f1YldawiugygYiIxJOTvkxQwN3Xufuj7j4QuAJYSugXePQzNSIiIlKjVGQI4/tdQxiLiIjUeBrCWEREJM5pCGMREZE4pyGMRURE4pyGMBYREYlzJV4mcPcFhMYHAMDM6hAaNTAN6IkuGYiIiNR4FRm18PkqrpuIiIjEQFn6GZgCtAV+DuQCTYFnzGxkVVZMREREYqMsTxOk8d2ohT0iRy0EZldp7URERKTKlSUZmMJ3LQh7zewRQqMWHi9+EREREakpSr1M4O5/iuGohSIiIhJjZe10CAB3/xY46VELRURE5NRTpoGKRERE5PSlZEBERCTOKRkQERGJc0oGRERE4pySARERkTgX82TAzAaa2Xoz22hmE4uINzGzuWa2ysw+NLOUYP5FZrYy4nUwGFERM3vQzLZHxAbFeLdERERqrHI9WlhRZlYb+A3Qn1DXxh+Z2Tx3XxtRbBKw0t2HmlliUP5Kd18PpEesZzswN2K5R919agx2Q0RE5LQS65aBHsBGd/886LPgJWBIoTJJBCMluvs6IMHMWhQqcyWwyd23VHWFRURETnexTgZaA9sipnODeZE+Bm4AMLMehAZJalOozC3Ai4Xm3RVcWphpZk0qr8oiIiKnt1gnA1bEPC80/TDQxMxWAncDK4D88ArMzgCuA16JWGYa0IHQZYSdwP8UuXGzcWaWbWbZe/bsOcldEBEROb3E9J4BQi0BF0RMtwF2RBZw94PAGAAzM2Bz8CpwDbDc3XdFLBP+bGa/A/5U1MbdfQYwAyAjI6NwEiIiIhKXYt0y8BHQyczaBb/wbwHmRRYws3OCGMBYYHGQIBQYQaFLBGbWKmJyKLC60msuIiJymoppy4C755vZXcBbQG1gpruvMbM7gvh0oDPwnJkdA9YCPyxY3swaEHoS4UeFVv3fZpZO6JJDThFxERERKYa5x2dreUZGhmdnZ1d3NURERGLCzJa5e0ZRMfVAKCIiEueUDIiIiMQ5JQMiIiJxTsmAiIhInFMyICIiEueUDIiIiMQ5JQMiIiJxTsmAiIhInFMyICIiEueUDIiIiMQ5JQMiIiJxTsmAiIhInFMyICIiEueUDIiIiMQ5JQMiIiJxTsmAiIhInFMyICIiEueUDIiIiMS5mCcDZjbQzNab2UYzm1hEvImZzTWzVWb2oZmlRMRyzOwTM1tpZtkR8881s7+a2YbgvUms9kdERKSmi2kyYGa1gd8A1wBJwAgzSypUbBKw0t1TgVuBxwvFL3f3dHfPiJg3EVjg7p2ABcG0iIiIlEGsWwZ6ABvd/XN3/xZ4CRhSqEwSoS903H0dkGBmLUpZ7xDg2eDzs8D1lVZjERGR01ysk4HWwLaI6dxgXqSPgRsAzKwH0BZoE8QceNvMlpnZuIhlWrj7ToDg/bwqqLuIiMhpqU6Mt2dFzPNC0w8Dj5vZSuATYAWQH8T6uPsOMzsP+KuZrXP3xWXeeCiBKEgi8sxsfblqX7JmwN5KXF880jGsHDqOFadjWHE6hhVX2cewbXGBWCcDucAFEdNtgB2RBdz9IDAGwMwM2By8cPcdwftuM5tL6LLDYmCXmbVy951m1grYXdTG3X0GMKNS9yhgZtmF7mOQctIxrBw6jhWnY1hxOoYVF8tjGOvLBB8BncysnZmdAdwCzIssYGbnBDGAscBidz9oZg3NrHFQpiFwNbA6KDcPGB18Hg28UcX7ISIictqIacuAu+eb2V3AW0BtYKa7rzGzO4L4dKAz8JyZHQPWAj8MFm8BzA01FlAHeMHd/xLEHgbmmNkPga3A8Fjtk4iISE0X68sEuPt8YH6hedMjPr8HdCpiuc+BtGLWuQ+4snJrWm5VcvkhzugYVg4dx4rTMaw4HcOKi9kxNPfC9++JiIhIPFF3xCIiInFOyUAlKK2LZSmZmV1gZgvN7FMzW2Nm91R3nWoqM6ttZivM7E/VXZeaKLiB+VUzWxecj72ru041jZn9W/B3vNrMXjSz+tVdp5rAzGaa2W4zWx0xL2Zd7SsZqKAydrEsJcsHfuLunYFewHgdw5N2D/BpdVeiBnsc+Iu7JxK6R0nHshzMrDXwYyDD3VMI3Sh+S/XWqsaYBQwsNC9mXe0rGai4snSxLCVw953uvjz4fIjQf8CFe6aUUphZG+Ba4OnqrktNZGZnAX2B3wO4+7fuvr9aK1Uz1QHONLM6QAMK9SUjRQs60Puy0OyYdbWvZKDiytLFspSRmSUAXYEPqrkqNdFjwE+B49Vcj5qqPbAHeCa41PJ00KeJlJG7bwemEnrEeydwwN3frt5a1Wgx62pfyUDFlaWLZSkDM2sE/AG4N+iJUsrIzAYDu919WXXXpQarA3QDprl7V+BrNAJquQTXtIcA7YDzgYZm9oPqrZWUhZKBiiu1i2UpnZnVJZQIzHb316q7PjVQH+A6M8shdKnqCjP7v+qtUo2TC+S6e0Gr1KuEkgMpu6uAze6+x92PAq8Bl1RznWqyXUEX+5TU1X5lUDJQcaV2sSwlC8ag+D3wqbs/Ut31qYnc/efu3sbdEwidg++4u36RlYO7fwFsM7OLgllXEuoFVcpuK9DLzBoEf9dXopswKyJmXe3HvAfC001xXSxXc7Vqmj7AKOCTYLRKgElBb5UisXQ3MDtI7D8nGDRNysbdPzCzV4HlhJ4SWoF6IiwTM3sR6Ac0M7NcYDIx7GpfPRCKiIjEOV0mEBERiXNKBkREROKckgEREZE4p2RAREQkzikZEBERiXNKBkREROKckgEREZE4p2RARKqEmZ1pZu+aWW0zS4gcp72c6znDzBYHo+CJSBVQMiAiVeU24DV3P1aRlQRDgy8Abq6UWolIFCUDIlJuZrbQzPoHn39pZk8UUWwkRfSlbmbtgyGCLw5aDNYFwwWvNrPZZnaVmf3dzDaYWY9gsdeD9YlIFVAyICInYzLwgJmNBLoC/xYZDPr2b+/uOYXmX0RodMox7v5RMLsj8DiQCiQC/wJcCtwPTArKrAYurpI9ERENVCQi5efui4NR6e4D+hVxKaAZsL/QvOaEWgpuLDSY12Z3/wTAzNYAC9zdzewTICHY3jEz+9bMGrv7ocrfI5H4ppYBESk3M+sCtAL+WcyX8zdA/ULzDgDbCI1SGemfEZ+PR0wf58QfLPWAIydbZxEpnpIBESkXM2sFzAaGAF+b2YDCZdz9K6C2mUUmBN8C1wO3mtm/lHObTYE97n70pCsuIsVSMiAiZWZmDYDXgJ+4+6fAQ8CDxRR/m9C1/zB3/xoYDPybmQ0px6YvB+aXu8IiUibm7tVdBxE5DZlZV+A+dx9VCet6Dfi5u6+veM1EpDC1DIhIlXD3FcBCM6tdkfUETya8rkRApOqoZUBERCTOqWVAREQkzikZEBERiXNKBkREROKckgEREZE4p2RAREQkzikZEBERiXNKBkREROLc/wNEctqVyZ53uQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M1 = timml.Model3D(kaq=K, z=[0., -b])\n",
    "timml.ImpLineDoublet1D(M1, xld=0.)\n",
    "timml.HeadLineSink1D(M1, xls=L, hls=BC3h, res=(BC3c+BC3c*dpar)/100., wh=1.)\n",
    "timml.StripAreaSink(M1, 0., L, R)\n",
    "M1.solve(silent=True)\n",
    "H1 = M1.headalongline(X, 0.).flatten()\n",
    "dhdBC3c = (H1-H0)/(BC3h*dpar)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., dhdBC3c, 'k--', mfc='none', label='perturbation sensitivity, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial h_{\\Gamma_1}$')\n",
    "plt.legend()\n",
    "plt.ylim(0.95, 1.05);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000000\n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f'% dhdBC3c[ocol])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2c. Numerical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAACdCAYAAACkc0HoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAlgElEQVR4nO3deVxV1d748c+XSUQERBwBAS1FBUVE0xweFZzNAZ+0zBxuWT5lt6zM7HavZb9b1GO3rl7TsknTrpHlVGblVNrNCSVF0XAAxFQEE82Jaf3+AM/jiKDA5sD3/XqdF+fsYe3v2R7P96y1115LjDEopZRSyv44WB2AUkoppW6NJnGllFLKTmkSV0oppeyUJnGllFLKTmkSV0oppeyUJnGllFLKTjlZHUBJ+fj4mMDAQKvDUEoppcpFXFxchjGmzvXW2V0SDwwMZNu2bVaHoZRSSpULEUm50TptTldKKaXslCZxpZRSyk7ZXXN6aXvmmWfYvXs3Dg4OtkfTpk2ZPn06AM899xypqak4ODjg6OiIg4MDLVq0YPLkyQDExMSQlZWFi4uL7dG0aVMGDRoEwBdffEFeXh4uLi5Uq1YNNzc3fH19ueOOOwDIzMykevXqVK9eHRGx5iQopZSyS1U+iZ89e5ZTp06Rn59ve3h4eNjWJyUlkZiYSF5enm19Tk6Obf2CBQtISkoiOzvbtmzw4MG2JD5+/HgyMjKuOObIkSP55JNPAPDz8+PChQsAuLm5UaNGDR5++GFeffVV8vLy6NGjB+7u7nh6euLh4YGnpydRUVH07NmTnJwcVq5ciaenp+3h4eGBl5cXTk5V/p9WKaUqvSr/TT9nzpwi1y9ZsqTI9QkJCQAYY8jNzSU7O5vLJ5XZvHkz58+fJzs7mwsXLnD+/Hl8fHxs6//3f/+Xc+fOcfbsWdvf0NBQAHJychARjh07xq+//kpWVhZZWVlUq1aNnj17kpmZyeDBg6+JKSYmhsmTJ5OWlka/fv2oXbu27eHt7c3QoUOJiIjgzJkz7Ny5E29vb9s6Tf5KKWU/9Bu7lIgIzs7OODs7X7G8cePGRe43YcKEG65zdXVl/fr11yzPz88HwNvbm23btpGVlcXp06dtSb5Tp04A5OXl0bhxYzIzM9mzZw+ZmZmcPHmS5s2bExERwa5du+jcufMVZXt5efHJJ58wYMAA9uzZw+zZs6lXr94Vj5YtW+Lu7l6c06KUUqoMaRK3Qw4OBf0RXVxcaNu27Q23CwgIYOnSpVcsM8bYfgS0aNGCVatWcfLkSTIzM8nMzCQjI4OAgAAAUlJSWLBgAadOnbqijB9//JEuXbqwePFipkyZQt26dWnYsCG+vr74+voyZswY6tSpw9mzZ3F0dMTV1bX03rxSSikbsbf5xCMiIozeJ16+Ll68SHp6OsePH+f48eN06tQJLy8v1q9fz7vvvsuxY8c4evQoaWlpnD17loMHDxIUFMQbb7zB5MmT8fb2xs/Pz5bk33zzTTw8PDh06BDnzp0jICBAa/ZKKXUDIhJnjIm47jpN4qq0GGM4ffo0NWvWxMHBgc2bN7N69WqOHDnCkSNHSEtL47fffiMlJQUXFxf+/Oc/M3PmTABq165NQEAAQUFBxMbG4uDgQEJCArm5uQQGBuLl5WXtm1NKKYtoElcV0t69e4mPjyc5OZmUlBRSUlI4d+6crR/AkCFDbJcDPDw8CAgIoG3btnz00UcA7NixAzc3N4KCgnBxcbHoXSilVNnSJK7s0u7du0lMTLQl+OTkZDw8PFiwYAEAbdu2Zfv27Tg4OODv788dd9xBZGQkU6ZMAeDAgQPUr1+fGjVqWPk2lFLqthSVxLVjm6qwWrZsScuWLW+4fs6cOezdu5f9+/ezf/9+Dhw4wMGDB23r7777btLT02nQoAFNmjShadOm9O/fn+joaKCg976jo2OZvw+llCormsSV3WrXrh3t2rW77jpjDDNnzrQl+P3797Ny5Uq8vLyIjo7m/Pnz1KpVi6CgIIKDg22Pzp0706RJk3J+J0opdWs0iatKSUQYNmzYNcvz8vIAyM7O5umnn2bv3r3s3buXr776itzcXP7xj38wceJEUlJSGD16NMHBwbRs2ZKQkBBCQkKoU+e6swEqpZQlNImrKuVS87mnpyevvvqqbXlOTg6HDh3C09MTgNOnT5OTk0NsbCy///67bbsvvviC6OhoDhw4wNq1awkJCaFly5ZXDNWrlFLlRZO4UoCzszNNmza1vQ4NDeWnn37CGMPx48fZtWsXCQkJtub7NWvW8Oijj9q2b9SoESEhIcyePZtGjRpx+vRpXF1dtde8UqpMae90pW5Bfn4+KSkpJCQk2B67du1iw4YNeHp68re//Y2YmBhatGhBWFgYbdq0ISwsjM6dO2tnOqVUiegtZkqVs40bN7Jy5Up27NhBfHw8x44dw93dnaysLBwcHJg9ezbp6emEhYURFhZGo0aNdCpapdR16S1mSpWzzp07XzG5zLFjx0hOTraNe79mzRq+/PJL24x3tWrVYuDAgXz88ccAHD9+nLp162piV0oVSZO4UuWgfv361K9f3/Z68eLFnD17lp07dxIfH098fPwVU9SGh4dz7tw5IiIiaNeuHREREXTo0IGGDRtaEb5SqoLS5nSlKpj8/Hzef/99tm3bxtatW21jyD/xxBPMmDGDnJwcXn31VSIiIoiIiKBevXpWh6yUKkN6TVwpO3b+/Hl27tyJp6cnwcHBJCYm0rJlS1tTvL+/P+3bt+eZZ56hY8eOFkerlCptek1cKTtWvXp17rrrLtvr5s2bc/r0aXbs2MHWrVvZunUrmzZt4syZMwCsXbuWKVOm0LFjRzp27EiHDh2045xSlVS5JnERcQS2AUeMMQNExBv4DAgEkoFhxpjfb1yCUgrA3d2dLl260KVLF9uySzXz/Px8XF1dee+99/jnP/8JQIMGDfjpp58ICgoiIyODGjVqUL16dUtiV0qVHodyPt6TQOJlr58H1hhj7gTWFL5WSt2CSzXtqKgofvjhB7KysoiLi2PWrFn07t0bf39/AKZNm4aHhwcRERE89dRTfPHFF6Snp1sZulLqFpXbNXER8QPmAX8Hni6sie8DuhljjopIA2C9MaZZUeXoNXGlbs+GDRv45ptv+M9//sPmzZu5cOECDRs2JC0tDRFh48aN+Pr6EhgYqE3wSlUAFeWa+NvAc0DNy5bVM8YcBShM5HXLMR6lqqTLm+EvXrzI9u3bSU9PtyXskSNHkpKSQsOGDW3bRkZGEhwcbGXYSqnrKJfmdBEZAKQbY+Jucf9HRGSbiGw7ceJEKUenVNVVrVo1OnbsyKBBg4CC6+orVqxg1qxZdO3alY0bNzJhwgTefvttoGAWuOnTp7NlyxZyc3MtjFwpBeXUnC4irwEPArmAK+ABfAm0Q5vTlaqwjDGkpKRgjCEoKIjExERatGgBQM2aNenatSvdu3fn3nvvpVGjRhZHq1TlVFRzernUxI0xU4wxfsaYQOA+YK0xZiSwHBhduNloYFl5xKOUKh4RITAwkKCgIKDg9rZjx46xaNEiRowYQVJSEs8++ywJCQkA7N69mxkzZrBr1y7y8/OtDF2pKsHq+8RjgFgReQhIBe61OB6l1E3Uq1eP4cOHM3z4cACOHDmCt7c3AN999x1PP/00AD4+PnTr1o3u3bszZswY3NzcLItZqcpKR2xTSpWqlJQU1q1bZ3tkZGTw+++/U61aNWJjY7l48SI9e/a8Yix5pdSN6bCrSilLGGM4duwYDRo0ACAyMpK1a9cCEBoaSq9evRgwYADdunWzMEqlKjbLr4krpaomEbElcIDvv/+euLg4XnvtNerUqcPMmTOZOXOmbf27777Ljh079Hq6UsWkNXGllGXOnTvHyZMn8fPzu6LGXqdOHaKioujVqxd9+/bVmdpUlaY1caVUheTm5oafnx9QMOf6kSNHmDdvHr1792bt2rWMHTuWVatWAXD8+HE2bNig96crdRmtiSulKiRjDLt27aJRo0Z4eXkxa9YsJkyYgJeXl62G3qdPH+0gpyq9ijLsqlJKFZuI0KpVK9vrkSNHUr9+fVauXMk333xDbGwsjo6OnDhxglq1apGenk7t2rVxdHS0MGqlypcmcaWUXfD09GTo0KEMHToUYwy//PILcXFx1KpVC4AHH3yQuLg4evfuTb9+/ejduzc+Pj4WR61U2dLmdKVUpfDll1+ybNkyvvnmG06cOIGI8OijjzJ79myrQ1PqtpRqc7qI1AAuGGPybjsypZQqJdHR0URHR5Ofn8/27dtZsWIFjRs3BuDMmTN06tSJ3r17M3DgQDp27IiTkzZEKvt305q4iDhQMN75AxRMWHIRqAacAFYC7xljkso4ThutiSulSurgwYM89thjrF27lpycHGrXrk3//v15/vnnad68udXhKVWk273FbB3QBJgC1DfG+Btj6gJdgE1AjIiMLLVolVKqlDVu3JhVq1aRkZHB559/Tr9+/fjqq69st6v99NNPzJ49m7S0NIsjVapkilMTdzbG5NzuNqVFa+JKqdKQm5uLo6MjIsKzzz7Lm2++CUB4eLitaV5r6aoi0LHTlVKqCMYY9u7dy7Jly1i2bBmbNm0iICCAQ4cOISIcOXKEhg0bIiJWh6qqoFJJ4iISAfwFCKCgQ5wAxhjTqsgdS5kmcaVUWTty5AiHDh2ic+fO5OXl0aBBA9zd3W019A4dOuDgoANeqvJRWsOuLgQ+AoYC9wADCv8qpVSl4uvrS+fOnQHIy8sjJiaG5s2bM2PGDDp16oSfnx+LFi2yOEqlSpbETxhjlhtjDhljUi49yiwypZSqAFxcXPjTn/7E119/zYkTJ/j000+5++67bcO9xsXF2dZnZ2dbHK2qakrSnB4J3A+soeA2MwCMMV+WTWjXp83pSqmKZNGiRYwfP56srCy8vLwYMmQIw4YNIyoqSu9FV6WitK6JLwCCgd3Apcl+jTHmT6USZTFpEldKVTTZ2dmsXr2azz77jKVLl2KMIT09HVdXV5KTk/Hz89OErm5ZaSXxXcaY0FKN7BZoEldKVWQXL15k9+7dhIeHA9CyZUtOnDjB0KFDGTZsGF27dtVJWlSJlFbHtk0i0qKUYlJKqUqpWrVqtgRujOGVV16he/fuzJ8/nx49euDr68vcuXMtjlJVFiVJ4p2BeBHZJyI7RWSXiOwsq8CUUsreiQjR0dF89tlnpKenExsbS+fOnfHy8gIgNTWVp59+mm3btmFvY3aoiqEkzekB11te3j3UtTldKVVZLFmyhPvuu4/s7GyaNWvGiBEjGDFiBHfccYfVoakK5LauiYuImJtsVJxtSosmcaVUZfL777+zePFiFi5cyA8//ICjoyPp6el4e3uTl5en18/V7U+AIiJPiEijqwp1EZEeIjIPGF0agSqlVFVTq1Ytxo0bx/r160lNTeXTTz/F29sbgH79+tG3b18++eQTzpw5Y3GkqiIqThLvA+QB/xaR30Rkj4gcBJIouG/8LWPMx2UYo1JKVQn+/v4MGzYMKOgU1759exITExk1ahT16tXj/vvvZ8OGDRZHqSqSmyZxY8wFY8w7xphOFIybHgmEG2MCjDHjjDHxZR2kUkpVNSLCK6+8wsGDB9mwYQOjR4/mu+++Y+vWrQCcPXuWLVu2aIe4Kk5nMVNKKTuRnZ1Nbm4ubm5uzJ8/n9GjRxMcHMzo0aN58MEH8fX1tTpEVQZK6z5xpZRSFnJxccHNzQ2AQYMGMXfuXHx8fJgyZQqNGjWid+/e/PHHHxZHqcqTJnGllLJDnp6ePPzww2zYsIGkpCT+8pe/4Orqiru7OwAffvghGzdu1Ob2Sq4k94n3AB4ATgEJwE4gwRhzsaj9Sps2pyulVNFycnLw9/fn+PHjNGnShFGjRjFq1CgCAwOtDk3dgtJqTl8AfAVsAhoDf6NgMpTiBOAvIutEJFFEdovIk4XLvUXkexFJKvxbqwTxKKWUug5nZ2f279/PvHnzaNSoEVOnTiUoKIh//etfVoemSllJauI/GmO63tJBRBoADYwx20WkJhAHDAbGACeNMTEi8jxQyxgzuaiytCaulFIlk5KSwieffEJ0dDQtWrTgxx9/5NNPP+Whhx4iIiICEbE6RFWE26qJi8h8EXkK+FlEnrmVAIwxR40x2wufnwESAV9gEDCvcLN5FCR2pZRSpSggIIAXX3yRFi0K5rBKTExk/vz5tG/fnrCwMGbMmMHJkyctjlLdiuI0p88r3K4+8KCIpIjIchF5RUTuLekBRSQQaANsBuoZY45CQaIH6pa0PKWUUiXz6KOPcvToUebMmYOLiwtPPvkkbdu21U5wdqjE94mLiBPQAmgNtDLGTCrBvu7AD8DfjTFfisgpY4zXZet/N8Zcc11cRB4BHgFo1KhR25SUcp1zRSmlKrVffvmF1NRU7rnnHvLy8ujevTs9e/Zk7Nix+Pn5WR1elXdbE6CUYhDOFHSM+9YY84/CZfuAbsaYo4XXzdcbY5oVVY5eE1dKqbJz/PhxHnjgAdasWYODgwN9+vThoYce4p577sHZ2dnq8Kokywd7kYJeEx8AiZcSeKHl/N/kKaOBZeURj1JKqeurV68eq1ev5sCBA0yZMoX4+HiGDh3K6tWrAbTJvYIpl5q4iHQGNgC7gPzCxS9QcF08FmgEpAL3GmOK7F2hNXGllCo/ubm5rF69mp49e+Lo6MgLL7zAjh07GD9+PP3798fJycnqECu9omri5XL2jTEbgRvdwxBZHjEopZQqOScnJ/r06WN7XbduXXbu3MngwYPx9fVl3LhxPPzwwzpuu0V02FWllFLF9tRTT5GSksKSJUsICQnhpZdeYtKk/+vfrM3t5UuTuFJKqRJxcnJi8ODBrFq1iv379/Pyyy8DsHv3bpo0aUJMTAzp6ekWR1k1aBJXSil1y5o0acKdd94JwIULFwgICGDKlCn4+flx3333sX79eq2dlyFN4koppUpF27ZtWbduHYmJiTz++ON8++239O/fn9OnTwPa1F4WNIkrpZQqVcHBwbz11lv89ttvfP/993h6emKMITIykokTJ5KUlGR1iJWGJnGllFJlonr16tx9990AnD9/ngYNGjBr1iyaNm1K3759+frrr8nPz79JKaoomsSVUkqVOTc3NxYuXEhqaiovv/wyv/zyCwMGDGDRokVWh2bXNIkrpZQqN/Xr1+dvf/sbKSkpfPbZZwwZMgSAOXPmMG7cOHbu3GlxhPZFk7hSSqly5+zszLBhw6hevTpQMGb7woULad26NV27diU2Npbc3FyLo6z4NIkrpZSy3NSpU0lLS2P69OmkpaUxfPhwHnjgAavDqvA0iSullKoQvL29eeaZZ0hKSmL58uU88cQTABw9epQJEybw66+/WhxhxaNJXCmlVIXi6OjIPffcQ+fOnQH4+eefmTt3Ls2aNWPAgAGsWbNG7zkvpElcKaVUhRYdHU1qaipTp05ly5YtREVF0aZNGy5evGh1aJYrl6lIS9P1piLNyckhLS2NCxcuWBSVUsXn6uqKn58fzs7OVoeilN25cOEC//73v9m3bx8xMTEALFiwgJ49e1KvXj2LoysbRU1FWimS+KFDh6hZsya1a9dG5EYzniplPWMMmZmZnDlzhqCgIKvDUcru/fbbb/j7++Pk5MT999/PxIkTad26tdVhlaqiknilaE6/cOGCJnBlF0SE2rVra6uRUqWkYcOG7N27l3HjxvH5558TFhZGZGRklekEVymSOKAJXNkN/awqVbruvPNO/vWvf5GWlsbrr7/O4cOH8fb2BiA1NZXs7GyLIyw7lSaJW23GjBk0b968QtzX+NJLLzF9+vTbKmP9+vUMGDCgyG1OnTrFO++8Y3v922+/8d///d+3ddyK6tL4zyVVGv8WSqniqVWrFs899xz79u3Dx8cHYwzDhg2jcePGTJ8+3TabWmWiSbyUvPPOO6xcuZKFCxcWa/vSGonIGGPZBAJXJ/GGDRuyePFiS2IpK3l5eQD85z//sTgSpVRxXd7aNXXqVJo2bcqkSZPw9/fnueee48iRIxZGV7o0iZeC8ePHc/DgQQYOHMhbb73FyZMnGTx4MK1ataJDhw62sYBfeuklHnnkEXr16sWoUaPo16+fbV2bNm2YNm0aAH/96195//33+eOPP4iMjCQ8PJzQ0FCWLVsGQHJyMs2bN+exxx4jPDycw4cP8/e//51mzZoRFRXFvn37rhvn559/TkhIiG1YQyjoTzB27FhCQ0Np06YN69atu2a/q2uTISEhJCcn8/zzz3PgwAHCwsKYNGkSycnJhISEFFnuxx9/THR0NH369OHOO+/kueeeu+n57datG5MnT6Z9+/Y0bdqUDRs22MqaMGGCbbsBAwawfv16ANzd3Zk8eTJt27YlKiqKLVu20K1bNxo3bszy5cuBggQ9adIk2rVrR6tWrXj33XeBglaI7t27M2LECEJDQ23lXfLGG28QGhpK69atef755wGYO3cu7dq1o3Xr1gwdOpRz587d9H0ppcqWiNC3b1/Wrl3L1q1b6du3L2+++abtO8DeOnZfj5PVAZSFbt26XbNs2LBhPPbYY5w7d45+/fpds37MmDGMGTOGjIyMa5qELyWGG5kzZw6rVq1i3bp1+Pj48MQTT9CmTRuWLl3K2rVrGTVqFPHx8QDExcWxceNGqlevTkxMDBs2bCAwMBAnJyd++uknADZu3MjIkSNxdXVlyZIleHh4kJGRQYcOHRg4cCAA+/bt46OPPuKdd94hLi6ORYsWsWPHDnJzcwkPD6dt27bXxDlt2jS+/fZbfH19OXXqFACzZs0CYNeuXezdu5devXoVu0NITEwMCQkJtveWnJxsW1dUufHx8ezYsYNq1arRrFkznnjiCfz9/Ys8Vm5uLlu2bGHlypW8/PLLrF69usjtz549S7du3Xj99dcZMmQIL774It9//z179uxh9OjRDBw4kA8++ABPT0+2bt3KxYsX6dSpE7169QJgy5YtJCQkXNOD/JtvvmHp0qVs3rwZNzc3Tp48CRTcxzpu3DgAXnzxRT744APbaFNKKetFRESwaNEiXnvtNerXrw/Au+++y/Lly5k0aRLdunWzy/4qWhMvAxs3buTBBx8EoEePHmRmZpKVlQXAwIEDbQP+d+nShR9//JGNGzfSv39//vjjD86dO0dycjLNmjXDGMMLL7xAq1atiIqK4siRIxw/fhyAgIAAOnToAMCGDRsYMmQIbm5ueHh42BL91Tp16sSYMWOYO3eurZn48liDg4MJCAgolV6dRZUbGRmJp6cnrq6utGjRgpSUlJuWFx0dDUDbtm2v+LFwIy4uLvTp0weA0NBQ/uu//gtnZ2dCQ0Nt+3/33XfMnz+fsLAw7rrrLjIzM0lKSgKgffv2170FbPXq1YwdOxY3NzcAW+eZhIQEunTpQmhoKAsXLmT37t03jVEpVf6CgoJs38EODg7ExcXRo0cP2rVrR2xsrO270V5Uypp4UTVnNze3Itf7+PjctOZ9M9drorn0C69GjRq2Ze3atWPbtm00btyYnj17kpGRwdy5c2216IULF3LixAni4uJwdnYmMDDQdmvS5eVcXn5R5syZw+bNm/n6668JCwsjPj6+WM1JTk5OV1x3L87tUUWVW61aNdtzR0fHYvUPuLTP5dsXFZezs7PtnDg4ONj2d3BwsO1vjGHmzJn07t37imOtX7/+mvN7+fu63rkeM2YMS5cupXXr1nz88ce3/RlSSpW9Rx55hFGjRjF//nzefPNNhg8fzrBhw/jss8+sDq3YtCZeBrp27Wrr4LZ+/Xp8fHzw8PC4ZjsXFxf8/f2JjY2lQ4cOdOnShenTp9OlSxcAsrKyqFu3Ls7Ozqxbt+6GNdauXbuyZMkSzp8/z5kzZ1ixYsV1tztw4AB33XUX06ZNw8fHh8OHD18R66+//kpqairNmjW7Yr/AwEC2b98OwPbt2zl06BAANWvW5MyZMzc9Bzcq92qjRo1iy5YtRW5zdVzx8fHk5+dz+PDhEu0L0Lt3b2bPnk1OTo4tzrNnzxa5T69evfjwww9t17wvNaefOXOGBg0akJOTU+zOjUop67m6uvLII4+wZ88eFi9ezP/8z/8AkJGRwdtvv80ff/xhcYRFq5Q1cau99NJLjB07llatWuHm5sa8efNuuG2XLl1Ys2YNbm5udOnShbS0NFsSf+CBB7jnnnuIiIggLCyM4ODg65YRHh7O8OHDCQsLIyAgwLb/1SZNmkRSUhLGGCIjI2ndujXBwcGMHz+e0NBQnJyc+Pjjj6+oKQMMHTrU1uzcrl07mjZtCkDt2rXp1KkTISEh9O3bl8cff9y2z2OPPXbTcq+2c+dOGjRoUOQ2l+vUqRNBQUGEhoYSEhJCeHh4sfcFePjhh0lOTiY8PBxjDHXq1GHp0qVF7tOnTx/i4+OJiIjAxcWFfv368eqrr/LKK69w1113ERAQQGho6A1/3CilKiZHR0eGDh1qe7106VImTpzIK6+8wpNPPsmECRNsl88qkkox7GpiYiLNmze3KCJVGk6fPs1DDz3E559/bnUo5UI/s0pVfD///DOvvfYaK1aswN3dnfHjx/P666/j4FC+jdiVfthVZf88PDyqTAJXStmHjh07snz5cnbu3MnAgQM5ePCgLYGfOHHC4ugKaBJXSimlinDprpPY2FgAkpKS8PX1ZcSIEbaxPqyiSVwppZQqBkdHRwC8vLx46qmnWLFiBa1bt2bAgAFs2rTJkpgqTRK3t2v7qurSz6pS9q1OnTq88cYbpKamMm3aNDZt2kS3bt3IzMws91gqRRJ3dXUlMzNTvxxVhXdpPnFXV1erQ1FK3aZatWrx17/+leTkZL766itq165d7jFYfouZiPQB/gk4Au8bY2JKWoafnx9paWkVpqOBUkVxdXXFz8/P6jCUUqXE3d2dqKgoS45taRIXEUdgFtATSAO2ishyY8yekpTj7Ox83SEylVJKqcrM6ub09sB+Y8xBY0w2sAgYZHFMSimllF2wOon7Aocve51WuEwppZRSN2F1Er/erB3X9E4TkUdEZJuIbNPr3koppVQBqzu2pQGXTyTtB/x29UbGmPeA9wBE5ISI3HzuyuLzATJKsbyqSM9h6dDzePv0HN4+PYe3r7TPYcCNVlg6drqIOAG/ApHAEWArMMIYU26TMYvIthuNSauKR89h6dDzePv0HN4+PYe3rzzPoaU1cWNMrohMAL6l4BazD8szgSullFL2zOrmdIwxK4GVVsehlFJK2RurO7ZVBO9ZHUAloOewdOh5vH16Dm+fnsPbV27n0O7mE1dKKaVUAa2JK6WUUnaqSidxEekjIvtEZL+IPG91PPZGRPxFZJ2IJIrIbhF50uqY7JWIOIrIDhH5yupY7JGIeInIYhHZW/h57Gh1TPZGRCYW/j9OEJF/i4jO0lMMIvKhiKSLSMJly7xF5HsRSSr8W6usjl9lk/hl47b3BVoA94tIC2ujsju5wDPGmOZAB+BxPYe37Ekg0eog7Ng/gVXGmGCgNXouS0REfIE/AxHGmBAK7ha6z9qo7MbHQJ+rlj0PrDHG3AmsKXxdJqpsEkfHbb9txpijxpjthc/PUPDFqcPmlpCI+AH9gfetjsUeiYgH0BX4AMAYk22MOWVpUPbJCaheOH6HG9cZeEtdyxjzI3DyqsWDgHmFz+cBg8vq+FU5ieu47aVIRAKBNsBmi0OxR28DzwH5FsdhrxoDJ4CPCi9JvC8iNawOyp4YY44A04FU4CiQZYz5ztqo7Fo9Y8xRKKjsAHXL6kBVOYkXa9x2dXMi4g58ATxljDltdTz2REQGAOnGmDirY7FjTkA4MNsY0wY4Sxk2X1ZGhddsBwFBQEOghoiMtDYqVRxVOYkXa9x2VTQRcaYggS80xnxpdTx2qBMwUESSKbik00NEFlgbkt1JA9KMMZdagRZTkNRV8UUBh4wxJ4wxOcCXwN0Wx2TPjotIA4DCv+lldaCqnMS3AneKSJCIuFDQiWO5xTHZFRERCq5DJhpj/mF1PPbIGDPFGONnjAmk4DO41hijNaASMMYcAw6LSLPCRZHAHgtDskepQAcRcSv8fx2Jdg68HcuB0YXPRwPLyupAlg+7ahUdt71UdAIeBHaJSHzhshcKh9JVqjw9ASws/EF+EBhrcTx2xRizWUQWA9spuOtkBzpyW7GIyL+BboCPiKQBU4EYIFZEHqLgB9K9ZXZ8HbFNKaWUsk9VuTldKaWUsmuaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslOaxJVSVxCR6iLyQ+Ec54GXz5NcwnJcROTHwlmxlFJlQJO4UupqfwK+NMbk3U4hhVP8rgGGl0pUSqlraBJXqgoRkXUi0rPw+f8TkRnX2ewBrjPWs4g0Lpzqs11hDX1v4bSfCSKyUESiROQnEUkSkfaFuy0tLE8pVQY0iStVtUwF/iIiD1Aw//vEy1cWjj3e2BiTfNXyZhTMVjfWGLO1cPEdwD+BVkAwMALoDDwLvFC4TQLQrkzeiVKq6k6AolRVZIz5sXCWqqeBbtdpMvcBTl21rA4FNfOhV00SdMgYswtARHYDa4wxRkR2AYGFx8sTkWwRqWmMOVP670ipqk1r4kpVISISCjQALt4gqZ4HXK9algUcpmDWustdvOx5/mWv87myglANuHCrMSulbkyTuFJVhIg0ABYCg4CzItL76m2MMb8DjiJyeSLPBgYDo0RkRAmPWRs4YYzJueXAlVI3pElcqSpARNyAL4FnjDGJwCvASzfY/DsKrm3bGGPOAgOAiSIyqASH7g7o/PJKlRGdT1wpdQURaQM8bYx5sBTK+hKYYozZd/uRKaWupjVxpdQVjDE7gHUi4ng75RT2dF+qCVypsqM1caWUUspOaU1cKaWUslOaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslOaxJVSSik7pUlcKaWUslP/H/j0MwEZPo+rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import flopy\n",
    "\n",
    "dpar = 1e-3\n",
    "nrow, ncol = 1, int(L)\n",
    "M0 = flopy.modflow.Modflow(modelname='model', exe_name='../mf2005.exe')\n",
    "flopy.modflow.ModflowDis(M0, nlay=1, nrow=nrow, ncol=ncol, nper=1, delr=1., delc=1., top=0., botm=-b, steady=True, \n",
    "                         perlen=1., nstp=1)\n",
    "flopy.modflow.ModflowBas(M0, ibound=np.ones([nrow, ncol]), strt=BC3h*np.ones([nrow, ncol], dtype=float))\n",
    "flopy.modflow.ModflowLpf(M0, hk=K, vka=-999., ss=-999., sy=-999., ipakcb=53)\n",
    "flopy.modflow.ModflowRch(M0, nrchop=1, rech=R, ipakcb=53)\n",
    "flopy.modflow.ModflowGhb(M0, stress_period_data={0: [0, 0, ncol-1, BC3h, BC3c]})\n",
    "flopy.modflow.ModflowPcg(M0, hclose=1e-6, rclose=1e-6)\n",
    "flopy.modflow.ModflowOc(M0, stress_period_data={(0,0): ['save head', 'save budget']})\n",
    "M0.write_input()\n",
    "success, buff = M0.run_model(silent=True)\n",
    "H0 = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "hBC3 = H0[-1]\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., H0, 'k--', mfc='none', label='forward solution, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel('$h$ (m)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAAChCAYAAABJR/xvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhS0lEQVR4nO3deXxU9b3/8dcHCKuigBuCElQsS0jCIqsSrCiIWJYfuFxEglLAoj/1ior6aO293F71QtW6XBSXqogLIgpVFJWCYksLQSllEQFBCFJkUcQFJPC5f8xJnGQmk0nIJJPwfj4e5zFzvud7vudzvplkPjnL95i7IyIiIlKcGpUdgIiIiCQ3JQsiIiISk5IFERERiUnJgoiIiMSkZEFERERiUrIgIiIiMSVNsmBmjc3sXTNbH7w2KqbejWa2ysxWm9lNYeW/NbNtZrYimPqHLbvDzDaY2Toz61sBuyMiIlJtJE2yAEwEFrh7K2BBMF+ImaUBvwS6ABnAADNrFVblAXfPDKZ5wTptgSuAdkA/4H/NrGZid0VERKT6SKZkYSDwbPD+WWBQlDptgL+5+/funge8DwyOo92X3P2Au28CNhBKNkRERCQOyZQsnOzu2wGC15Oi1FkF9DKzJmZWH+gPnBa2/HozW2lmT4edxmgGbA2rkxuUiYiISBxqVeTGzOw94JQoi+6KZ313X2tm9wHvAt8C/wDygsVTgUmAB6+/B64BLFpTxcQ3BhgD0KBBg06tW7eOJywREZFqYfny5bvc/cSi5RWaLLh7n+KWmdkOM2vq7tvNrCnwZTFtPAU8Fazz34SOFODuO8LaegJ4I5jNpfDRh+bAF8W0PQ2YBtC5c2fPycmJc89ERESqPjP7PFp5Mp2GmAuMDN6PBOZEq2RmJwWvpwNDgBeD+aZh1QYTOmWR3+4VZlbHzFoCrYCl5R69iIhINVWhRxZKcC8w08yuBbYAwwDM7FTgSXfPvxXyVTNrAhwExrv7V0H5/5hZJqFTDJuBsQDuvtrMZgJrCJ2yGO/uhypml0RERKo+0yOqo9NpCBEROdqY2XJ371y0PJmOLIiIxO3gwYPk5uayf//+yg5FpMqpW7cuzZs3JyUlJa76ShZEpErKzc3l2GOPJTU1FbNoNz2JSDTuzu7du8nNzaVly5ZxrZNMFziKiMRt//79NGnSRImCSCmZGU2aNCnVUTklCyJSZSlRECmb0v7uKFkQEakkr7/+OmvWrCn1eosWLWLAgAGlWufBBx/k+++/L5jv378/X3/9dam3XZEee+wxnnvuOQCeeeYZvvjipyFyRo8eXWLf9ejRA4DNmzfzwgsvJC7QChDeF6V1zDHHHPH2lSyIiFSCvLy8MiULeXl5JVeKomiyMG/ePI4//vgytVVRxo0bx9VXXw1EJgtPPvkkbdu2jbn+X//6V6DqJwt5eXmF+qIyKFkQESmDzZs307p1a0aOHEl6ejpDhw4t+DJevnw5WVlZdOrUib59+7J9+3YAevfuzZ133klWVhb33Xcfc+fO5dZbbyUzM5ONGzfSu3dv8m/Z3rVrF6mpqUDoi3LYsGFceumlXHTRRQB88803DB48mLZt2zJu3DgOHz4MwHXXXUfnzp1p164dd999NwAPPfQQX3zxBeeffz7nn38+AKmpqezatQuA+++/n7S0NNLS0njwwQcL9q9Nmzb88pe/pF27dlx00UX88MMPEf3wyiuvkJaWRkZGBr169QLg0KFD3HrrrZxzzjmkp6fz+OOPA6EjIr1792bo0KG0bt2a4cOHk3/7/sSJE2nbti3p6elMmDABgN/+9rdMmTKFWbNmkZOTw/Dhw8nMzOSHH34o6KupU6dy2223FcTzzDPPcMMNNwA//Uc9ceJEFi9eTGZmJg888ADnnXceK1asKFinZ8+erFy5stifday4w/sxJyeH3r17F8Q+cuRILrroIlJTU5k9eza33XYb7du3p1+/fhw8eDDuz8of/vCHgr4A2LBhA3369CEjI4OOHTuyceNGvv32Wy644AI6duxI+/btmTMn6riGZefumqJMnTp1chFJXmvWrKnU7W/atMkB//DDD93dfdSoUT558mT/8ccfvXv37v7ll1+6u/tLL73ko0aNcnf3rKwsv+666wraGDlypL/yyisF81lZWb5s2TJ3d9+5c6e3aNHC3d3/+Mc/erNmzXz37t3u7r5w4UKvU6eOb9y40fPy8rxPnz4F7eTXycvL86ysLP/HP/7h7u4tWrTwnTt3Fmwrfz4nJ8fT0tL822+/9X379nnbtm39o48+8k2bNnnNmjX9448/dnf3YcOG+fTp0yP6IS0tzXNzc93d/auvvnJ398cff9wnTZrk7u779+/3Tp06+WeffeYLFy70hg0b+tatW/3QoUPerVs3X7x4se/evdvPPvtsP3z4cKF27r77bp88eXJE34TPf/nll37mmWcWlPfr188XL17s7u4NGjQo6K9LLrmkoM4zzzzjN954o7u7r1u3zkv6e19c3EX7ddmyZZ6VlVUQe8+ePf3HH3/0FStWeL169XzevHnu7j5o0CB/7bXXSvVZCe+LLl26+OzZs93d/YcffvDvvvvODx486Hv37nX30GfnzDPPLOjP/H4oKtrvEJDjUb4TdeukiFQL+f/Rhbvsssv41a9+xffff0///v0jlmdnZ5Odnc2uXbsYOnRooWWLFi0qcZunnXYaPXv2BOCqq67ioYceol+/fqxatYoLL7wQCP2X3bTpT6PRX3755aXYq59ceOGFNG7cuGC+S5cunHHGGQBceeWVfPjhhwwdOpSZM2cybdo08vLy2L59O2vWrCE9Pb3Ydj/88EMGDx5MgwYNABgyZAiLFy/mF7/4BS1btiQzMxOATp06sXnz5oj1e/bsSXZ2NpdddhlDhgwB4J133mHlypXMmjULgL1797J+/Xpq165Nly5daN68OQCZmZls3ryZbt26UbduXUaPHs0ll1xSqusxTjzxRM444wz+9re/0apVK9atW1fwMynOsGHDmDRpEpMnT+bpp58mOzu7xO1Ei/vcc8+Nuc7FF19MSkoK7du359ChQ/Tr1w+A9u3bs3nzZtatW1fqz8q+ffvYtm0bgwcPBkLjJUBo3JE777yTDz74gBo1arBt2zZ27NjBKadEe3Zj6SlZEBEpo6JXlJsZ7k67du1YsmRJ1HXyv5SjqVWrVsHphKK3tRVdL9q2N23axJQpU1i2bBmNGjUiOzu7xNvjPMYovnXq1Cl4X7NmzainIR577DH+/ve/8+abb5KZmcmKFStwdx5++GH69u1bqO6iRYsi2szLy6NWrVosXbqUBQsW8NJLL/HII4/w5z//OWbc4S6//HJmzpxJ69atGTx4cIlX+tevX58LL7yQOXPmMHPmTOIZrTda3BD7Z5a/To0aNUhJSSmIq0aNGuTl5ZXps1Lcz2vGjBns3LmT5cuXk5KSQmpqarkOWKZkQUSqhVhHAurXrx9z+QknnBDXkYSitmzZwpIlS+jevTsvvvgi5557Lj/72c/YuXNnQfnBgwf59NNPadeuXcT6xx57LPv27SuYT01NZfny5XTp0qXgv/LiLF26lE2bNtGiRQtefvllxowZwzfffEODBg047rjj2LFjB2+99VbBEZf8bZ1wwgmF2unVqxfZ2dlMnDgRd+e1115j+vTpcffBxo0b6dq1K127duVPf/oTW7dupW/fvkydOpWf//znpKSk8Omnn9KsWbNi2/j2228Ljv5069aNs846q8S+CjdkyBB+97vf0aJFC+6777641h09ejSXXnop5513XsERm6VLl/LII4+U6q6D/J/ZxRdfzKuvvhr3ekCpPiv5GjZsSPPmzXn99dcZNGgQBw4c4NChQ+zdu5eTTjqJlJQUFi5cyOefR314ZJnpAkcRkTJq06YNzz77LOnp6ezZs4frrruO2rVrM2vWLG6//XYyMjLIzMwsuCq/qCuuuILJkyfToUMHNm7cyIQJE5g6dSo9evQouGiuON27d2fixImkpaXRsmVLBg8eTEZGBh06dKBdu3Zcc801hQ7HjxkzhosvvrjgAsd8HTt2JDs7my5dutC1a1dGjx5Nhw4d4u6DW2+9lfbt25OWlkavXr3IyMhg9OjRtG3blo4dO5KWlsbYsWNj3sWxb98+BgwYQHp6OllZWTzwwAMRdbKzsxk3blzBBY7hGjVqRNu2bfn888/p0qVLxLrp6enUqlWLjIyMgrY7depEw4YNGTVqVEG9LVu2UK9evbj3HeDuu+/mxhtv5LzzzqNmzZqlWrc0n5Vw06dP56GHHiI9PZ0ePXrwr3/9i+HDh5OTk0Pnzp2ZMWMGrVu3LlUsJdGDpIqhB0mJJLe1a9fSpk2bStv+5s2bGTBgAKtWraq0GKTsvvjiC3r37s0nn3xCjRqh/5tvvfVWRowYEfMaj+ok2u9QcQ+S0pEFERE5qjz33HN07dqV3/3udwWJAsDkyZOPmkShtHTNgohIGaSmpuqoQhV19dVXV+oAR1WRjiyIiIhITEoWRKTK0jVXImVT2t8dJQsiUiXVrVuX3bt3K2EQKSV3Z/fu3QUDOsUjaa5ZMLPGwMtAKrAZuMzdv4pS70bgl4ABT7j7g0H5y8DPgmrHA1+7e6aZpQJrgXXBsr+5+7hE7YeIVIzmzZuTm5vLzp07KzsUkSqnbt26BSNSxiNpkgVgIrDA3e81s4nB/O3hFcwsjVCi0AX4EXjbzN509/XufnlYvd8De8NW3ejumYneARGpOCkpKbRs2bKywxA5KiTTaYiBwLPB+2eBQVHqtCF0ZOB7d88D3gcGh1ew0HialwEvJi5UERGRo0cyJQsnu/t2gOD1pCh1VgG9zKyJmdUH+gOnFalzHrDD3deHlbU0s4/N7H0zOy8RwYuIiFRXFXoawszeA6I9AuuueNZ397Vmdh/wLvAt8A+g6BiiV1L4qMJ24HR3321mnYDXzaydu38TJb4xwBiA008/PZ6QREREqr0KTRbcvU9xy8xsh5k1dfftZtYU+LKYNp4CngrW+W8gN6yNWsAQoFNY/QPAgeD9cjPbCJwNRIzl7O7TgGkQGu651DsoIiJSDSXTaYi5wMjg/UhgTrRKZnZS8Ho6ocQg/ChCH+ATdw9PIE40s5rB+zOAVsBn5R69iIhINZVMd0PcC8w0s2uBLcAwADM7FXjS3fsH9V41sybAQWB8kdsrryDywsZewH+aWR5wCBjn7nsSuB8iIiLVip46WQw9dVJERI42euqkiIiIlImSBREREYlJyYKIiIjEpGRBREREYlKyICIiIjEpWRAREZGYyjzOgpldHbxtCBxw9yfKJyQRERFJJkdyZOFkoLW7PwKcWk7xiIiISJI5khEc33f3pcH7j8ojGBEREUk+pTqyYGYT8t+HJQq4+5/KMygRERFJHqU9DXFGQqIQERGRpFXa0xAtzKx/eIG7zyvHeERERCTJlDZZmAmcGDavp1CJiIhUczGTBTP7NfC9u/8ewN2frZCoREREJGmUdGRhBJBZtNDMRgMnuvs9iQhKREREkkdJFzj+4O7fRymfDlyVgHhEREQkyZSYLJhZ06KF7n4AyEtMSNXXoUOHKmSd8lTR2z+S7VV2X1WmZNr38oglmfZHREpOFn4PzDGzFuGFZnYScDhhUVVDZkatWrX45JNPMLNC029+85uIstmzZxesU3RZ+LR3796Isjlz5sRcJ3zKj63oNHTo0Ijtz507N6428//QFy1/+OGHI8qaNWtWaB/ytzdx4sS49yFarCNGjIio17Fjx4iyyy+/PKJs+fLlEWX3339/1L7btm1bRHnXrl255JJLIspzcnLo0aNHobKbbrqJsWPHxr2vRacRI0bQsmXLiM/JlClTIupee+21EWVr166N2u7jjz8etf82bNgQUf7oo4+yfv36iJ/hV199xTnnnFOo7po1a5g/f37B/BNPPBF1+/lt/PWvf42rH2bMmBG1PDs7O6Ls3nvvjSh74YUXIspWrlwZtc3bb7+90Hy9evW4/vrrI+r1798/6mdGpEpy95gTMBLYBbwB/Bfw38B6YHhJ61blqVOnTl6eCN05knTTkiVLyr3N7777rlT7fNVVV1V6P4RPvXv3jqvegQMHKj3WZJgmTZrkZ599dlx1p02bVmj+rLPOqvT4K3I6fPhwuf5dESlvQI5H+U4scVAmD90B0ZLQbZMpwH7gSnefUdK6pWFmw8xstZkdNrPOMer1M7N1ZrbBzCaGlTc2s3fNbH3w2ihs2R1B/XVm1rc8467qunfvXu5t1q9fv1T1n3/++XKP4UgsWrQornpr1qxJbCBVRM2aNdm/f39cdRcvXlxo/rjjjktESEmrRo0aZT6KVNap6NGs6jxJ4sQ1gqO773P359z9dnf/T3fPSUAsq4AhwAfFVTCzmsCjwMVAW+BKM2sbLJ4ILHD3VsCCYJ5g+RVAO6Af8L9BOyJSDp588km2bNkSV90aNQr/ybn00ksTEZKEWbJkSWWHUGEqO1mJNs2fP59TTjkFMyM9PR0z44ILLqBx48aYGa1bt6ZJkya0a9eu0HoNGzbEzDjmmGM488wzI9p94okKftBztMMN4RNQB/g1MA94BbgDOLOk9co6AYuAzsUs6w7MD5u/A7gjeL8OaBq8bwqsK1onmJ8PdC8pjqPlNET41K9fv7jr1qlTp9hldevW9ZNPPrnS9ydR0z333OOA165du8S6V199ddTyW265JWp5WlpaRNntt99e6fscbUpNTT3iNnJzcyPKOnbsGLXuqaeeGleb5557rg8cOLBc97V3nKemipsaNWpU6T8vTdVvuv76633Xrl3l/V0V9TREPF/eDwFPAhnAA8AUQv/9X1XSumWZiJ0sDAWeDJsfATwSvP+6SN2vgtdHwmMFngKGFtP+GCAHyDn99NPLs/99/PjxhX7II0aMiPgSeOONN7xHjx6Fyk466SQ/7bTTCuazsrJ87ty5UT840a4/2L59e0TZ0qVLC96Hf+Ht2bMnou7UqVMj/tBNmDDBt23b5iNGjChU/vLLL3u9evV87NixPnbs2ILy/C/X8Om9997znj17ekpKipuZA964cWNfv359oXo33HCDb926NWL9HTt2RJQ9/fTT3qFDh0JlU6ZM8aysrEJl99xzj48cObJQ2WOPPeYnnnhiobJu3bpFxJ7/ubj55psL7eMFF1zgq1evLlS3cePGvmfPHh80aFCh8lGjRvmmTZu8Ro0ahcrnzZvn27Zti9gv98hk8/nnn48omzlzZkTZ3LlzI75833nnHf/1r38dUXfevHl+zjnnOFCQ7N10001RY3r//ff9lVde8Ro1avhFF11U0HcTJkzwzMzMQn0wcOBAnz9/vgOF+v3++++Pum/RPrOLFi3ytWvXFiqrX7++b9iwIaLu22+/7c8991zEZ/Ott96KqPvqq69GlBX9bAwcOND37NkTkRyOHz/eZ8+eHbH+unXr4ior7VT081m0L4pbdtVVV3mbNm3i3k6tWrXiqlf081uWqXHjxhFlxxxzTFz7lZGRccTbr4jpgw8+KEh0O3XqVOzPtejfqfzpP/7jP/xXv/qVn3DCCX7LLbd4enq6N2vWzKdPn+5Lliwp1++p4PcxarJgoWXFM7P3gYHu/rWZ/cXde5pZHWCRu5fqhLeZvQecEmXRXe4+J6izCJjgUU51mNkwoK+7jw7mRwBd3P0GM/va3Y8Pq/uVuzcys0eBJe7+fFD+FDDP3V+NFWvnzp09JycRZ1tERESSk5ktd/eI6wbjeTbEZH66tmGXmd0PfEIZbp109z6lXaeIXOC0sPnmwBfB+x1m1tTdt1tobIgv41hHREREShDzAkcz6w686e57gqJhhA7THxO8r2jLgFZm1tLMahO6cHFusGwuods8CV7nhJVfYWZ1zKwl0ApYWoExi4iIVGklHVkYCTxqZp8CbwNvu/sLiQjEzAYDDxN6quWbZrbC3fua2amErlPo7+55ZnY9oYsUawJPu/vqoIl7gZlmdi2whSCZcffVZjYTWENo1Mnx7q7h4UREROJU4jULAGbWmtDtin2B44CFhJKHv1TXL15dsyAiIkeb4q5ZiHechU/c/QF37wf8HPiQ0H/ufy/fMEVERCTZxHOBI2b2c2A48DWhwZNWErpj4UDiQhMREZFkEFeyADwPjA/qpwODCI2IeFZiwhIREZFkEW+ysMHdXwvev5KoYERERCT5lHTr5HNmdhOwxMxuqZiQREREJJmUdIHjs0GdU4ARZva5mc01s0nBaIoiIiJSzcU8DeHuCwg9wREAM6tF6GmPGUBXdEpCRESk2ivxmoXgORC3EXri43fAR8BMd5+e4NhEREQkCcQzzsJkoAWhRz3nAk2AP5rZ8EQGJiIiIskhnrshMvjpqZNdwp86CcxIaHQiIiJS6Sr0qZMiIiJS9ZR4GsLd30iip06KiIhIBYt3UCYA3P1HICFPnRQREZHkFNeDpEREROTopWRBREREYlKyICIiIjEpWRAREZGYlCyIiIhITEmTLJjZMDNbbWaHzaxzjHr9zGydmW0ws4lh5ZPN7BMzW2lmr5nZ8UF5qpn9YGYrgumxCtgdERGRaiNpkgVgFTAE+KC4CmZWE3gUuJjQA62uNLO2weJ3gTR3Twc+JTQ8db6N7p4ZTOMSEr2IiEg1lTTJgruvdfd1JVTrAmxw98+CMR9eAgYG67/j7nlBvb8BzRMXrYiIyNEjaZKFODUDtobN5wZlRV0DvBU239LMPjaz983svEQGKCIiUt2UagTHI2Vm7wGnRFl0l7vPiaeJKGVeZBt3AXn89JCr7cDp7r7bzDoBr5tZO3f/Jkp8Y4AxAKeffnoc4YiIiFR/FZosuHufI2wiFzgtbL458EX+jJmNBAYAF7i7B9s8ABwI3i83s43A2YSecVE0vmnANIDOnTt70eUiIiJHo6p2GmIZ0MrMWppZbeAKYC6E7pIAbgd+4e7f569gZicGF0ZiZmcArYDPKjxyERGRKippkgUzG2xmuUB34E0zmx+Un2pm8wCCCxivB+YDa4GZ7r46aOIR4Fjg3SK3SPYCVprZP4BZwLiwp2iKiIhICSw4Wi9FdO7c2XNyIs5UiIiIVFtmttzdI8Y6SpojCyIiIpKclCyIiIhITEoWREREJCYlCyIiIhKTkgURERGJScmCiIiIxKRkQURERGJSsiAiIiIxKVkQERGRmJQsiIiISExKFkRERCQmJQsiIiISk5IFERERiUnJgoiIiMSkZEFERERiUrIgIiIiMSlZEBERkZiULIiIiEhMSZMsmNkwM1ttZofNrHOMev3MbJ2ZbTCziWHlvzWzbWa2Ipj6hy27I6i/zsz6JnpfREREqpNalR1AmFXAEODx4iqYWU3gUeBCIBdYZmZz3X1NUOUBd59SZJ22wBVAO+BU4D0zO9vdDyVgH0RERKqdpDmy4O5r3X1dCdW6ABvc/TN3/xF4CRhYwjoDgZfc/YC7bwI2BO2IiIhIHJImWYhTM2Br2HxuUJbvejNbaWZPm1mjONcRERGRGCr0NISZvQecEmXRXe4+J54mopR58DoVmBTMTwJ+D1xTwjpF4xsDjAlmvzWzko50lNYJwK5ybvNooz48curDI6c+PHLqw/JR3v3YIlphhSYL7t7nCJvIBU4Lm28OfBG0vSO/0MyeAN4oaZ0o8U0Dph1hjMUysxx3L/biTSmZ+vDIqQ+PnPrwyKkPy0dF9WNVOw2xDGhlZi3NrDahCxfnAphZ07B6gwldMEmw/Aozq2NmLYFWwNIKjFlERKRKS5pkwcwGm1ku0B1408zmB+Wnmtk8AHfPA64H5gNrgZnuvjpo4n/M7J9mthI4H7g5WGc1MBNYA7wNjNedECIiIvEz96in7yUBzGxMcKpDykh9eOTUh0dOfXjk1Iflo6L6UcmCiIiIxJQ0pyFEREQkOSlZqADFDVEt8TGz08xsoZmtDYYEv7GyY6qqzKymmX1sZm+UXFuiMbPjzWyWmX0SfCa7V3ZMVY2Z3Rz8Lq8ysxfNrG5lx5TsgvGDvjSzVWFljc3sXTNbH7w2itXGkVCykGBhQ1RfDLQFrgyGoJb45QG3uHsboBswXn1YZjcSujhYyu4PwNvu3hrIQP1ZKmbWDPj/QGd3TwNqErqzTWJ7BuhXpGwisMDdWwELgvmEULKQeGUZolrCuPt2d/8oeL+P0B9njcJZSmbWHLgEeLKyY6mqzKwh0At4CsDdf3T3rys1qKqpFlDPzGoB9Slm7Bv5ibt/AOwpUjwQeDZ4/ywwKFHbV7KQeBpuuhyZWSrQAfh7JYdSFT0I3AYcruQ4qrIzgJ3AH4PTOU+aWYPKDqoqcfdtwBRgC7Ad2Ovu71RuVFXWye6+HUL/VAEnJWpDShYSL+7hpiU2MzsGeBW4yd2/qex4qhIzGwB86e7LKzuWKq4W0BGY6u4dgO9I4KHf6ig4rz4QaEnoScANzOyqyo1KSqJkIfHiHm5aimdmKYQShRnuPruy46mCegK/MLPNhE6F/dzMnq/ckKqkXCDX3fOPbM0ilDxI/PoAm9x9p7sfBGYDPSo5pqpqR/7oxcHrl4nakJKFxCt2iGqJj5kZoXPEa939/sqOpypy9zvcvbm7pxL6DP7Z3fXfXCm5+7+ArWb2s6DoAkKjw0r8tgDdzKx+8Lt9AbpItKzmAiOD9yOBeB7IWCYV+iCpo5G755lZ/hDVNYGnw4aolvj0BEYA/zSzFUHZne4+r/JCkqPYDcCMIPn/DBhVyfFUKe7+dzObBXxE6E6nj0ngA/yqCzN7EegNnBA8GuFu4F5gppldSygJG5aw7WsERxEREYlFpyFEREQkJiULIiIiEpOSBREREYlJyYKIiIjEpGRBREREYlKyICIiIjEpWRAREZGYlCyISKUxs3pm9r6Z1TSzVDNbVcZ2apvZB8FTDEWknClZEJHKdA0w290PHUkjwePfFwCXl0tUIlKIkgURSQgzW2hmFwbv/8vMHopSbThRxrM3szOCR0CfExxx+CR4HPQqM5thZn3M7C9mtt7MugSrvR60JyLlTMmCiCTK3cBdZjYc6ADcHL4weLbCGe6+uUj5zwg9YXSUuy8Lis8C/gCkA62BfwPOBSYAdwZ1VgHnJGRPRI5yOr8nIgnh7h8ETxX8d6B3lFMNJwBfFyk7kdCRhv9X5IFrm9z9nwBmthpY4O5uZv8EUoPtHTKzH83sWHffV/57JHL00pEFEUkIM2sPNAUOFPPl/QNQt0jZXmAroSeNhjsQ9v5w2PxhCv/TUwfYX9aYRSQ6JQsiUu7MrCkwAxgIfGdmfYvWcfevgJpmFp4w/AgMAq42s38r5TabADvd/WCZAxeRqJQsiEi5MrP6wGzgFndfC0wCfltM9XcIXXtQwN2/AwYAN5vZwFJs+nxgXqkDFpESmbtXdgwicpQysw7Av7v7iHJoazZwh7uvO/LIRCScjiyISKVx94+BhWZW80jaCe6seF2Jgkhi6MiCiIiIxKQjCyIiIhKTkgURERGJScmCiIiIxKRkQURERGJSsiAiIiIxKVkQERGRmJQsiIiISEz/B4kwovznT/gfAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M1 = M0\n",
    "flopy.modflow.ModflowGhb(M1, stress_period_data={0: [0, 0, ncol-1, BC3h, BC3c+BC3c*dpar]})\n",
    "M1.write_input()\n",
    "success, buff = M1.run_model(silent=True)\n",
    "H1 = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "dhdBC3c = (H1-H0)/(BC3c*dpar)\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., dhdBC3c, 'k--', mfc='none', label='perturbation sensitivity, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial C_{\\Gamma_3}$')\n",
    "plt.legend()\n",
    "plt.ylim(-1.05, -0.95);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-0.999451\n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f'% dhdBC3c[ocol])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Adjoint sensitivity\n",
    "\n",
    "<br>\n",
    "    \n",
    "$\\begin{align}\n",
    "\\quad\\quad \\frac{\\partial h(x')}{\\partial h_{\\Gamma_3}}=\\oint\\limits_{\\Gamma_3}\\psi^*_1(x)\\ [h_{\\Gamma_3}-h(x)]\\; dx\\ =\\ \\psi^*_1(\\Gamma_3)\\ [h_{\\Gamma_3}-h(\\Gamma_3)]\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Governing equation:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad K\\ b\\ \\frac{d\\psi^*_1}{dx}+\\frac{1}{2\\ K\\ b}\\delta(x-x')=0\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Boundary conditions: \n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad \\psi^*_1(x)=0\\ ,\\quad\\quad\\quad\\quad x=0\\cr\n",
    "\\end{align}$\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad(\\rm{undefined\\ on\\ }\\Gamma_3)\\cr\n",
    "\\cr\\end{align}$\n",
    "\n",
    "Closed-form solution:\n",
    "\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3a. Analytical\n",
    "$\\begin{align}\n",
    "\\quad\\quad\\rm{Not\\ available}\n",
    "\\cr\\end{align}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b. Semi-analytical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAACdCAYAAACtmOGSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkaUlEQVR4nO3de1xVZdr/8c8FaoSdNM0sCURTU1DwkCimpFaeMm3SashDP82xzBwtDSvUX6b2POZYOs2gNamN2kHTKaun7ACajWZalAc8AipZSajoqCTC9fzBdj9ydKvAgs31fr14sde6117ry9rg5b322vctqooxxhhjKj8fpwMYY4wxpnRYUTfGGGO8hBV1Y4wxxktYUTfGGGO8hBV1Y4wxxktYUTfGGGO8RDWnA1yqOnXqaFBQkNMxjDHGmHKxefPm31S1blFtlb6oBwUFsWnTJqdjGGOMMeVCRPYV12aX340xxhgvYUXdGGOM8RKV/vJ7abv77rv5/fff863r06cPTzzxBLm5ufTs2bPQc+677z4eeeQRTpw4wR/+8IdC7YMGDSI6Opr09HQGDx5cqH3EiBH079+f/fv3M3LkyELtTzzxBD169GDXrl2MHTu2UPvTTz9N586d+eGHH3j22WcLtU+ePJl27dqxYcMGpk2bVqh9xowZhISEkJCQwF/+8pdC7bNnz6ZRo0Z8/PHHxMXFFWqfN28e9evXZ8WKFSxatKhQ+6JFi7jmmmtYunQp7777bqH2t99+Gz8/P9544w1WrVqVr01EWLFiBQCvvvoqX3zxRb52f39/Fi9eDMCsWbP497//na/92muvZf78+QBMnz6d7777Ll/7jTfeyCuvvALApEmTSEpKytfeuHFjZsyYAcCECRNITU3N1x4SEsKkSZOAvNfp119/zdfetm1bxo8fD8Cf/vQnMjMz87V36tSJxx9/HIAhQ4Zw+vTpfO3du3dn2LBh5Obm8tBDD1FQnz59+OMf/8jJkyd55JFHCrXfd9999O/fn8OHD/PEE08Uao+OjqZnz5789NNPxMTEFGofNmwYUVFRJCcnM2XKlHxtAQEBPP7449SvX7/Q84wxzrCiXsB//vOfQkX93OVjx44Ves7ZdlXl6NGjhdqzsrLc7YcPHy62PScnh99++63Y/WdnZ3Po0KFi20+fPs0vv/xSqP1socjKyuLgwYPFtp84cYK0tLRC7dnZ2R61Hzt2jP379xdqz8nJAeDo0aOFiiLknReAjIwMkpOT87WJiPtxeno6e/bsyddes2ZN9+Off/6ZXbt25Wu/7rrr3I8PHDjAjh078rWf+9qmpqayffv2fO0+Pv93MWvPnj2Fnu/v7+9+vHPnzkI/f506ddyPt27dSkZGRr72gIAA9+PExET378JZTZo0cT8u6t6RsLAwIO8cb9y4sVB7REQEkPdzbtiwoVB7165dATh16lSh/xBB3n8aAI4fP866devc61WVAwcO0KlTJyvqxlQgUtkndGnbtq3ajXLGlL/9+/cTEBCAiBAbG8vPP//MxIkTadSokdPRjPFqIrJZVdsW1WbvqRtjLspNN93kvpJy5swZFi9eTNOmTRk8eHChKxrGmPJhRd0Yc8lmzJhBSkoKY8aM4b333qN58+a8/PLLTscypsqxom6MKRX169dn1qxZpKamEhMTQ1RUFADJycls3rzZ2XDGVBFW1I0xpapu3bpMnz7dfRPfjBkzaNu2Lb1792b9+vXOhjPGy1lRN8aUqVmzZjF9+nQ2btxIx44d6d69O2vXrnU6ljFeyYq6MaZMXXXVVUycOJHU1FReeukltm7d6h57QFWp7J/AMaYisaJujCkXNWvW5MknnyQlJcU9kE18fDwRERGsWrXKirsxpcCKujGmXF1++eVcc801QN6ASOnp6fTt25fw8HCWL19Obm6uswGNqcSsqBtjHNOrVy927tzJwoULOXXqFAMGDOCOO+5wOpYxlZYVdWOMo6pXr86QIUPYvn07b731FkOHDgXyBrRZunSpexhiY8z5WVE3xlQIvr6+PPDAAwwaNAiAVatWER0dTZMmTZg3b16hORmMMYVZUTfGVEj9+vXjww8/pF69eowcOZJGjRoxZ84c67kbUwIr6saYCklE3APWrF69muDgYOLi4vD19QWwG+qMKYIVdWNMhSYi3HHHHaxdu5Z169bh4+PDsWPHaNasGdOmTSs0R70xVZkVdWNMpVG7dm0Ajh49SpMmTXjuuecIDAxk0qRJHD582OF0xjjPiroxptK56aab+PDDD9m8eTNdu3Zl6tSpBAYG8tNPPzkdzRhHOVbURcRXRL4XkQ9dy7VF5DMR2e36XsupbMaYyqF169asWLGCLVu2EBMTw4033gjAsmXLOHjwoMPpjCl/TvbUxwBJ5yzHAF+o6s3AF65lY4w5r5CQEJ599lkAMjMzGTp0KMHBwYwaNYp9+/Y5nM6Y8uNIUReRBkBv4PVzVt8DLHI9XgT0K+dYxhgvcPXVV7NlyxYGDx7Ma6+9RuPGjRk+fDhpaWlORzOmzDnVU38ZmACc+5mUeqr6M4Dr+3UO5DLGeIHg4GDmz5/P3r17GTlyJO+88w6nT58G8kaqM8ZblXtRF5E+wCFV3XwJ+xghIptEZFN6enoppjPGeJOAgADmzp3LwYMHCQ4OBuDee+/l/vvv58cff3Q4nTGlz4meeiTQV0RSgbeBriKyGPhVROoDuL4fKm4HqjpfVduqatu6deuWR2ZjTCV25ZVXAnkD1oSGhvI///M/tGrVin79+rFp0yaH0xlTesq9qKvqRFVtoKpBwAPAl6r6EPABMMS12RDg/fLOZozxbj4+PkybNo3U1FSmTJnCmjVraNeuHW+++abT0YwpFRXpc+ovAneIyG7gDteyMcaUutq1azN58mT27dvHiy++SJ8+fQBYu3YtCQkJqKrDCY25OFLZf3nbtm2rdvnMGFMaevTowaeffkqnTp147rnnuPPOOxERp2MZk4+IbFbVtkW1VaSeujHGOGrlypXMnTuX1NRUevToQUREBF9++aXTsYzx2EUVdRGpKSK+pR3GGGOcdPnll/P444+zZ88e5s+fT3p6OsnJyQBkZ2fbzHCmwvOoqIuIj4j8UUQ+EpFDwA7gZxHZJiIzReTmso1pjDHl57LLLuORRx5h586dDBmSd/9uXFwcoaGhLF261D7rbiosT3vq8UAjYCJwvaoGqOp1wG3ABuBFEXmojDIaY4wjqlevTvXq1QFo2LAhIkJ0dDS33HILCxYsIDs72+GExuTn0Y1yIlJdVUv87fVkm7JgN8oZY8pLbm4u77//PlOnTuX7779nwIABvPvuu07HMlXMJd8o50mxdqKgG2NMefLx8aF///5s3ryZjz76iDFjxgBw6NAh5syZw8mTJx1OaKq68xZ1EblDRF4TkTDX8ogyT2WMMRWYiNCrVy8iIyMBWL58OWPGjKFhw4bMnDmT//znPw4nNFWVJz31x4DxwEMi0hUIK9NExhhTyTz22GOsWbOGVq1aMWHCBAIDA5k+fboNYmPKnSdFPV1Vj6rqU8CdQLsyzmSMMZVO586dWb16NevXr6djx44kJia6B645ceKEw+lMVeFJUf/o7ANVjQFskGRjjClGREQEq1atYsmSJQAkJSVxww038PTTT3PoULHzVBlTKs5b1FX1/QLLc8sujjHGeIezH4W77LLL6NWrFzNnziQoKIixY8dy8OBBh9MZb3XJw8SKyNelEcQYY7xRcHAwb731FklJSQwcOJC5c+cSEhJid8qbMlEaY7/fUAr7MMYYr9a0aVMWLlzIrl27iIuLw9/fH1Vl+vTp7Nmzx+l4xkt4OkzsXBEZISIdROTKAs12e6cxxngoODiYgQMHArB3716ef/55mjZtyqBBg9ixY4fD6Uxl52lPfQvQkrw5zlNFJEVEPhCRaUDBIm+MMcYDjRs3JiUlhbFjx7JixQqaN2/O/fffzy+//OJ0NFNJXdR86iLSgLwiHwqEqqpj477bMLHGGG+Qnp7O7NmzWb58OYmJifj7+5OZmcnVV1/tdDRTwZQ0TKynY793ADZoBRxJwYq6Mcab5OTk4Ovry5kzZ2jWrBlNmjThueeeo2PHjk5HMxXEJY/9DgwBvhORt0VkqIhcX3rxjDHGnOXr6wvAmTNnGD58ON9++y2RkZF069aN+Ph4G6XOlMjTCV1Gqmo4MAWoBSwUkfUiMl1EOouIb1mGNMaYqsbPz4+YmBhSU1OZNWsW27dvp2vXrnz++edORzMVmKeX36up6pkC6y4Hbgd6Ah2KuxRQ1uzyuzGmKsjKyuKdd95h0KBB+Pj4sHDhQmrXrs3dd9/tHo7WVA2lcfl9Y8EVqnpKVT9W1dFOFXRjjKkq/Pz8GDJkCD4+Pqgqr776Kvfccw/h4eEsW7aM3NxcpyOaCsDTom7/DTTGmApCRFi/fj2LFi3i1KlTDBw4kJCQEOLj452OZhxWzcPt6orIuOIaVfUvpZTHGGOMB6pVq8bgwYOJjo5m2bJlTJs2DT8/PwCOHj1KzZo13ePPm6rD0566L3AFeQPNFPVljDHGAb6+vjzwwAP8+OOPdOjQAYCnn36am2++mbi4OH7//XeHE5ry5GlP/WdVfb5MkxhjjLlo594s179/f3744QceffRRXnjhBSZMmMDw4cPx9/d3MKEpD/aeujHGeJkePXqwfv16PvvsMxo1asSYMWN45plnnI5lyoGnPfWu59tARKQijjhnjDFVkYjQvXt3unfvztq1awkMDARg48aNrF69mtGjR9sQtF7I0576ShEZLSI3nbtSRGqISFcRWUTeqHPGGGMqmM6dO7uL+urVq4mNjSUwMJBJkyaRkZHhcDpTmjwt6j2AHOAtETkoIttFJBnYDTwIzFbVhWWU0RhjTCl57rnn+O677+jWrRtTp04lKCiI6dOnOx3LlBKPLr+rahbwN+BvIlIdqAOcUtWjZZjNGGNMGQgPD+e9995j27ZtTJs2zT1wTW5uLr/++iv169d3OKG5WJ721N1UNVtVf76Ugi4iASISLyJJIrJNRMa41tcWkc9EZLfre62LPYYxxpiStWjRgqVLl/Lss88CsHLlSoKCgnjsscfYt2+fw+nMxbjgol5KzgBPquotQAQwSkSaAzHAF6p6M/CFa9kYY0wZOvtxuDZt2jB06FBef/11GjduzLBhw9izZ4/D6cyFcKSou3r637keHweSgBuBe4BFrs0WAf2cyGeMMVVRUFAQ8+bNY+/evYwcOZIlS5bQu3dvm+61EnGqp+4mIkFAOPANUE9Vf4a8wg9c52A0Y4ypkgICApg7dy4pKSksWrQIEeHUqVMMHz6cH374wel4pgTnLeoiUtP1/YrSPrhrn+8Bf1bVYxfwvBEisklENqWnp5d2LGOMMUD9+vWJiIgAIDExkXfffZewsDDuuecebMrrismTnnotEXkc6FSaB3bdRf8esERVV7hW/yoi9V3t9YFDRT1XVeeraltVbVu3bt3SjGWMMaYIHTp0YN++fUyZMoW1a9fSrl07evbsyfHjx52OZs7hSVHvBgwFgkWkxMvhIvK0JweVvLsy/gEkFZjh7QP+bxCbIcD7nuzPGGNM2atVqxaTJ09m3759zJgxgxo1anDFFXkXcZOTk+299wpAzvciiMgtQHUgQFU/KtD27rmLQJjrzvXz7bMT8BWwBch1rX6GvPfV3wVuAvYDA1T1cEn7atu2rdplIGOMcc6hQ4cICgqidevWxMbGcuedd+abYMaULhHZrKpti2o7b09dVZNU9ceCBd3lmKoOdH0NAD73JJCqrlNVUdWWqhrm+vpYVTNUtZuq3uz6XmJBN8YY47yrrrqKmTNnsn//fnr06EH79u354IMPrOfugAu6+11E3hKRp0Wkp4jUBV4osMmzpRfNGGNMZeDn58eoUaPYs2cPr732GhkZGfTr14/du3c7Ha3KudCPtMUBJ4E/AJ8Bs0XkyrON1rM2xpiqq0aNGgwfPpydO3fy5Zdf0qRJEwBiYmJYvHgxZ86ccTih97vQon4TsAYYqaphwHLg/5d2KGOMMZVXtWrViIqKAiArK4tPPvmEQYMGccstt/DGG2+QnZ3tbEAvdqFFvQl5l9x3iMj3wJ3AvSLSzXU53hhjjHHz8/Pju+++Y+XKlVx11VUMGzaMm2++mfXr1zsdzStdaFGfo6p9VbUx0BVYDPgD0cDHpR3OGGNM5efj40O/fv3YtGkTH330EcHBwQQHBwOQkpLCyZMnHU7oPc77kbZ8G+f1zmuRN1b7TqAZsFNVx5RNvPOzj7QZY0zlpKpERkayd+9ennrqKR599FH3595N8S7pI23nUtVwoBHwFLAW+Cvw5CUnNMYYU+WICC+++CKtWrViwoQJBAYG8sILL5CZmel0tErrYuZTz1HVbaq6QlU/VFW7ndEYY8xF6dy5M6tXr2bDhg107NiR2NhY3nzzTadjVVrVnA5gjDHGtG/fnlWrVvH999/TtGlTABYvXsyWLVt48sknue46m7TTE45PvWqMMcacFR4ejr+/PwBbtmzhpZdeIigoiLFjx/LTTz85nK7is6JujDGmQvqv//ovkpKSuP/++5k7dy7BwcG89NJLTseq0KyoG2OMqbCaNGnCggUL2L17Nw8//DCNGzcG4OjRo+zZs8fhdBWPFXVjjDEVXsOGDYmLi6Nfv34AzJkzh6ZNm/LQQw+RlJTkbLgKxIq6McaYSmfEiBGMGzeOf/3rX7Ro0YKBAwfyww8/OB3LcVbUjTHGVDrXX389M2fOJDU1lYkTJ/LJJ5/wzDPPOB3LcRc0olxFVNSIctnZ2aSlpZGVleVQKmOc5+fnR4MGDahevbrTUYwpc0eOHCEzM5OgoCCSk5MZPXo0zzzzDJGRkU5HK3UljSjnlZ9TT0tL48orryQoKAgRcTqOMeVOVcnIyCAtLY2GDRs6HceYMlerVi1q1aoFwK5du9i4cSOdOnXi9ttvJzY2lqioqCpRD7zy8ntWVhbXXnttlXgBjSmKiHDttdfa1SpTJfXo0YPU1FRmzZpFUlISXbt2pVu3buTk5Dgdrcx5ZVEHrKCbKs/+BkxVVrNmTcaNG0dKSgp//etfiYyMxNfXF4B169ZR2d96Lo7XFvXKYOHChTz++OMAxMXFnXe840mTJvH555+XuE1CQgL//ve/LzlbamoqISEh591u+vTp+ZY7dux4yccuD56c70t17utb0jYHDx50Lw8fPpzt27df8LESEhLo06fPBT/PGG/n5+fHqFGjmDp1KgAbN27ktttuIzw8nGXLlpGbm+twwtJlRb2CGDlyJIMHDy5xm+eff57u3buXuE1pFXVPFSzq5XnsS+HJ+S4PBYv666+/TvPmzR1MZIx3a926NYsWLSIrK4uBAwcSEhLCkiVLOHPGO+Yms6JeRvr160ebNm1o0aIF8+fPd69fsGABTZo0oUuXLnz99dfu9VOmTHEPf5iYmEhERAQtW7akf//+HDlyBIChQ4eyfPlyAIKCgpg8eTKtW7cmNDSUHTt2kJqaSlxcHLNnzyYsLIyvvvoqX6Y1a9YQFhZGWFgY4eHhHD9+HFVl/PjxhISEEBoayjvvvFPoZynY4+zTpw8JCQnExMRw6tQpwsLCiI6OBnDPhVzcfhMSEoiKiuK+++6jWbNmREdHn/cy2LJlywgJCaFVq1Z07twZgJycHMaPH0+7du1o2bIl8+bNc++/S5cuDBw4kCZNmhATE8OSJUu49dZbCQ0NZe/evYXOt6ev3RVXXMGzzz5Lq1atiIiI4NdffwVg1apVtG/fnvDwcLp37+5ef9bx48dp2LAh2dnZABw7doygoCCWLVvGpk2biI6OJiwsjFOnThEVFcXZT3N88skntG7dmlatWtGtWzcgr5fRsWNHwsPD6dixIzt37izx3Blj8qtWrRqDBw9m27ZtvP322/j6+jJ69GhOnjzpdLTSoaqV+qtNmzZa0Pbt2/Mtd+nSpdDXq6++qqqqJ06cKLJ9wYIFqqqanp5eqM0TGRkZqqp68uRJbdGihf7222968OBBDQgI0EOHDunvv/+uHTt21FGjRqmq6uTJk3XmzJmqqhoaGqoJCQmqqhobG6tjxoxRVdUhQ4bosmXLVFU1MDBQ58yZo6qqr776qg4bNqzQfgrq06ePrlu3TlVVjx8/rtnZ2bp8+XLt3r27njlzRn/55RcNCAjQgwcPakpKirZo0UJVVRcsWODOqarau3dvjY+PV1XVmjVr5jvG2eXi9hsfH69XXXWVHjhwQHNycjQiIkK/+uqrEs9lSEiIpqWlqarqkSNHVFV13rx5OnXqVFVVzcrK0jZt2mhycrLGx8fr1VdfrQcPHtSsrCy94YYbdNKkSaqq+vLLL7vPZUnnqajXTlUV0A8++EBVVcePH+8+/uHDhzU3N1dVVV977TUdN25cofM2dOhQXblypTv72W26dOmi3377rfvYZ5cPHTqkDRo00OTk5HyZMjMzNTs7W1VVP/vsM7333ntVVTU+Pl579+5d6Gcp+LdgjMkvJydHd+7c6X7cq1cv/dvf/qZZWVkOJysesEmLqYnWUy8jc+bMcffoDhw4wO7du/nmm2+Iioqibt261KhRg/vvv7/Q8zIzMzl69ChdunQBYMiQIaxdu7bIY9x7770AtGnThtTU1PNmioyMZNy4ccyZM4ejR49SrVo11q1bx4MPPoivry/16tWjS5cufPvttxf/g7uUtN9bb72VBg0a4OPjQ1hY2HmzR0ZGMnToUF577TX33aurV6/mzTffJCwsjPbt25ORkcHu3bsBaNeuHfXr1+eyyy6jUaNG3HnnnQCEhoZ6dJ6Keu0AatSo4X7f+txznpaWxl133UVoaCgzZ85k27ZthfY5fPhwFixYAORdrXn44YdLzLBhwwY6d+7s/jha7dq1gbzfjwEDBhASEsLYsWOLPJYxxnM+Pj40adIEgPT0dI4cOcJjjz1GcHAwr7zySqXrwXvl59QLSkhIKLbN39+/xPY6deqU2F7c8T7//HPWr1+Pv78/UVFR7o8WleYdyZdddhkAvr6+Hr0fFBMTQ+/evfn444+JiIjg888/9+gO0GrVquW7mcSTj0mVtN+zucGz7HFxcXzzzTd89NFHhIWFkZiYiKoyd+5c7rrrrnzbJiQk5Nu/j4+Pe9nHx6fQsQ4cOMDdd98N5L3P3qxZs2Jfu+rVq7tfv3Nzjx49mnHjxtG3b18SEhKYMmVKoZ8hMjKS1NRU1qxZQ05OznlvQlTVIn9XYmNjuf3221m5ciWpqalERUWVuB9jjOfq1avH119/zZdffsnUqVP585//zPTp04mPj68097pYT70MZGZmUqtWLfz9/dmxYwcbNmwAoH379iQkJJCRkUF2djbLli0r9Nyrr76aWrVqud8P/+c//+nutXviyiuv5Pjx40W27d27l9DQUJ5++mnatm3Ljh076Ny5M++88w45OTmkp6ezdu1abr311nzPCwoKIjExkdzcXA4cOMDGjRvdbdWrV3e/V3wuT/Zb0MSJE1m5cmWRudu3b8/zzz9PnTp1OHDgAHfddRd///vf3cfetWsXJ06cOO/5KSggIIDExEQSExMZOXJksa9dSTIzM7nxxhsBWLRoUbHbDR48mAcffDBfL72416tDhw6sWbOGlJQUAA4fPlzoWAsXLvT45zTGeEZE6NatGwkJCaxdu5ZevXq5e/Lr16/n6NGjzgY8DyvqZaBHjx6cOXOGli1bEhsbS0REBAD169dnypQpdOjQge7du9O6det8zzvbM1u0aBHjx4+nZcuWJCYmMmnSJI+Pfffdd7Ny5coib5R7+eWX3TecXX755fTs2ZP+/fvTsmVLWrVqRdeuXfnv//5vrr/++nzPi4yMpGHDhoSGhvLUU0/lyz1ixAhatmzpvlHuLE/2W9CWLVuK3Gb8+PGEhoYSEhJC586dadWqFcOHD6d58+a0bt2akJAQ/vSnP5XK3avFvXYlmTJlCgMGDOC2226jTp06xW4XHR3NkSNHePDBB93rhg4dysiRI903yp1Vt25d5s+fz7333kurVq3cb9VMmDCBiRMnEhkZWSUG0jDGSbfddhsLFiygWrVqnD59mv79+xMYGEhsbCwZGRlOxyuSV479npSUxC233OJQooszevRoWrdufd73Wr3ZXXfdxaeffup0jDKzfPly3n//ff75z3+W2zEr49+CMRVVYmIiL7zwAu+99x41a9bkscce48knn6RevXrlmqOksd+tp14BxMbG8s0339C3b1+nozjKmwv66NGjiYmJITY21ukoxpiLFBYWxvLly9m6dSt9+/Zl1qxZbN261elY+VhP3RgvZn8LxpSd1NRUAgMDERFiYmLIzMwkJiaGwMDAMj2u9dSNMcaYUnbuTKC///47//jHP2jcuDHDhg1jz549jmSqcEVdRHqIyE4R2SMiMRe7n8p+BcKYS2V/A8aUn9mzZ5OcnMyjjz7K0qVLadq0KXFxceWeo0IVdRHxBV4FegLNgQdF5II/HOjn50dGRob9o2aqLHXNp+7n5+d0FGOqjAYNGjBnzhxSUlIYN26ce1jr8lTRBp+5FdijqskAIvI2cA9wQdNWNWjQgLS0NNLT08sgojGVg5+fHw0aNHA6hjFVzvXXX8/MmTMdOXZFK+o3AgfOWU4D2l/oTqpXr+4eXtMYY4ypKirU5XegqDFUC11DF5ERIrJJRDZZb9wYY4zJU9GKehoQcM5yA+BgwY1Udb6qtlXVtnXr1i23cMYYY0xFVtGK+rfAzSLSUERqAA8AHzicyRhjjKkUKtzgMyLSC3gZ8AXeUNVp59k+HdhXihHqAL+V4v6qIjuHl87OYemw83jp7BxeutI+h4GqWuRl6gpX1J0mIpuKG6nHeMbO4aWzc1g67DxeOjuHl648z2FFu/xujDHGmItkRd0YY4zxElbUC5vvdAAvYOfw0tk5LB12Hi+dncNLV27n0N5TN8YYY7yE9dSNMcYYL2FF3aW0ZoerykQkQETiRSRJRLaJyBinM1VWIuIrIt+LyIdOZ6mMROQaEVkuIjtcv48dnM5U2YjIWNff8VYReUtEbHYgD4jIGyJySES2nrOutoh8JiK7Xd9rldXxrahTerPDGc4AT6rqLUAEMMrO40UbAyQ5HaISewX4RFWbAa2wc3lBRORG4AmgraqGkDduyAPOpqo0FgI9CqyLAb5Q1ZuBL1zLZcKKeh737HCqeho4OzucuQCq+rOqfud6fJy8f0hvdDZV5SMiDYDewOtOZ6mMROQqoDPwDwBVPa2qRx0NVTlVAy4XkWqAP0UM2W0KU9W1wOECq+8BFrkeLwL6ldXxrajnKWp2OCtGl0BEgoBw4BuHo1RGLwMTgFyHc1RWwUA6sMD1FsbrIlLT6VCViar+BLwE7Ad+BjJVdbWzqSq1eqr6M+R1foDryupAVtTzeDQ7nPGMiFwBvAf8WVWPOZ2nMhGRPsAhVd3sdJZKrBrQGvi7qoYDJyjDy53eyPWe7z1AQ+AGoKaIPORsKuMJK+p5PJodzpyfiFQnr6AvUdUVTuephCKBviKSSt7bQF1FZLGzkSqdNCBNVc9eJVpOXpE3nusOpKhquqpmAyuAjg5nqsx+FZH6AK7vh8rqQFbU89jscKVARIS89zGTVPUvTuepjFR1oqo2UNUg8n4Pv1RV6yFdAFX9BTggIk1dq7oB2x2MVBntByJExN/1d90Nu9nwUnwADHE9HgK8X1YHqlZWO65MVPWMiDwOfMr/zQ63zeFYlVEkMAjYIiKJrnXPqOrHzkUyVdRoYInrP+nJwMMO56lUVPUbEVkOfEfep1q+x0aW84iIvAVEAXVEJA2YDLwIvCsiw8j7D9OAMju+jShnjDHGeAe7/G6MMcZ4CSvqxhhjjJewom6MMcZ4CSvqxhhjjJewom6MMcZ4CSvqxhhjjJewom6MMcZ4CSvqxpgSicjlIrLGNcd70LnzRF/gfmqIyFrXrF/GmDJgRd0Ycz7/D1ihqjmXshPXtMZfAPeXSipjTCFW1I2pwkQkXkTucD1+QUTmFLFZNEWMVS0iwa6pTdu5evA7XNOcbhWRJSLSXUS+FpHdInKr62n/cu3PGFMGrKgbU7VNBp4VkWggHBh7bqNr7PRgVU0tsL4pebPxPayq37pWNwZeAVoCzYA/Ap2Ap4BnXNtsBdqVyU9ijLEJXYypylR1rWsWrnFAVBGX2OsARwusq0tez/0PBSY+SlHVLQAisg34QlVVRLYAQa7j5YjIaRG5UlWPl/5PZEzVZj11Y6owEQkF6gO/F1NkTwF+BdZlAgfIm5XvXL+f8zj3nOVc8ncgLgOyLjazMaZ4VtSNqaJEpD6wBLgHOCEidxXcRlWPAL4icm5hPw30AwaLyB8v8JjXAumqmn3RwY0xxbKibkwVJCL+wArgSVVNAqYCU4rZfDV57427qeoJoA8wVkTuuYBD3w58fMGBjTEesfnUjTElEpFwYJyqDiqFfa0AJqrqzktPZowpyHrqxpgSqer3QLyI+F7Kflx30v/LCroxZcd66sYYY4yXsJ66McYY4yWsqBtjjDFewoq6McYY4yWsqBtjjDFewoq6McYY4yWsqBtjjDFewoq6McYY4yX+F3d4JNs9SlrXAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M0 = timml.Model3D(kaq=K, z=[0., -b])\n",
    "timml.ImpLineDoublet1D(M0, xld=0.)\n",
    "timml.HeadLineSink1D(M0, xls=L, hls=0., res=BC3c/100., wh=1.)\n",
    "timml.LineSink1D(M0, xls=float(ocol), sigls=-1.)\n",
    "M0.solve(silent=True)\n",
    "A = M0.headalongline(X, 0.).flatten()\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., A, 'k--', mfc='none', label='adjoint solution, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\psi^*_1$ (T/$L^2$)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAACdCAYAAAC+T6FmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAgiklEQVR4nO3de3xU9Z3/8deHcBPwx11KQQjZokJIIBgUZFdwQS4FQawXBF1AxcUFbKGCuGrh0daWLT68RLFULeCF6oLFigrKYg24RV0u0uW+grBLxMpFoHIP8Pn9MZPpZDIJk2Qyk5D38/HIg3PO53znfOZkwmfO5Xu+5u6IiIjIha9GshMQERGRxFDRFxERqSZU9EVERKoJFX0REZFqQkVfRESkmlDRFxERqSZqJjuBitasWTNPTU1NdhoiIiIJsW7dugPu3jxaLKFF38zmAoOBfe7eKUrcgKeB7wPHgdHuvj4YGxCMpQAvuvvMWLaZmprK2rVr4/QOREREKjcz+9/iYok+vT8fGFBCfCDQPvhzL/BrADNLAWYH4x2B282sY4VmKiIicoFJaNF391XANyWsMhR42QM+ARqZWUvgKmCHu3/h7qeB14PrioiISIwq2zX9VsCesPm84LJoy69OYF4iUkozZsxg9erVhZZ997vfZf78+QBMnTqVDRs2FIqnpaUxZ84cACZOnMj27dsLxTt16sQTTzwBwN13382ePXsKxbt168Zjjz0GwIgRIzhw4ECheK9evXj44YcBGDZsGMeOHSsUHzBgAJMnTw5Nnzt3rlB82LBh3HfffZw8eZIhQ4YUec8jRoxg9OjRHDp0iNtuu61I/J577uHWW2/lyy+/ZMyYMUXiEydO5IYbbmDnzp3cd999ReIPPvggffr0YePGjfz4xz8uEp8+fTo9e/bkv/7rv3jkkUeKxP/t3/6NrKwscnNz+cUvflEknpOTwxVXXMGyZct48skni8RffPFF2rRpw+LFi0O/p3ALFiygefPmLFiwgJdeeqlIfPHixTRo0IAXXniBRYsWFYkvW7aMlJQUcnJyeOeddwrF6tSpw9tvvw3AzJkz+eMf/1go3qhRIxYuXAhUrc9e586dmTVrVpF9UVEqW9G3KMu8hOXRX8TsXgKXB2jTpk18MhORUjlx4gRHjx4ttCy8yEaLnzhxIjR9/PjxUsdPnjxZaFuljZ86dSo0ffTo0SJFPzIe6fTp0wC4e5ni+fn5AJw9e7ZM8bNnzwJw5syZMsUL3m9+fn6Z4gVjuZw+fbpM8QKnTp0qEi9478XFa9WqFZquSp+98OWJYIkecMfMUoF3irmR7zdArru/FpzfDvQGUoEZ7t4/uPwhAHf/5fm2l52d7bqRT0REqgszW+fu2dFila2f/hLgnyygO3DE3b8C1gDtzaydmdUGhgfXFZFKqHv37tx+++3JTkNEIiS6y95rBI7cm5lZHjAdqAXg7nOApQS66+0g0GVvTDB2xswmAO8T6LI31903JzJ3EYnd8ePHC50KF5HKIaFF391L/OrvgWsN44uJLSXwpUBEKrlEXzYUkdhUttP7InKBCDxrS0QqExV9ERGRaqKyddkTkQvAzTffTIsWLZKdhohEUNEXkbibPn16slMQkSh0el9E4s7ddTOfSCWkoi8icZeRkcEtt9yS7DREJIKKvoiISDWhoi8icadT+yKVk4q+iFQI9dMXqXxU9EVERKoJddkTkbgbPXq0+umLVEIq+iISd1OmTEl2CiIShU7vi0jcHT16lBMnTiQ7DRGJoKIvInHXrVs3Ro8enew0RCSCir6IiEg1oaIvInGnfvoilVPCi76ZDTCz7Wa2w8ymRYlPMbMNwZ9NZnbWzJoEY7vNbGMwtjbRuYtI7NRPX6TySejd+2aWAswGrgfygDVmtsTdtxSs4+6zgFnB9W8AJrn7N2Evc527H0hg2iIiIheERHfZuwrY4e5fAJjZ68BQYEsx698OvJag3EQkTiZMmMB3vvOdZKchIhESXfRbAXvC5vOAq6OtaGb1gAHAhLDFDiw3Mwd+4+7PV1SiIlJ2EyZMOP9KIpJwiS760S7yFXfHzw3AnyJO7fd0971mdgnwH2a2zd1XFdmI2b3AvQBt2rQpb84iUkpff/01tWrVokmTJslORUTCJPpGvjzg0rD51sDeYtYdTsSpfXffG/x3H/AmgcsFRbj78+6e7e7ZzZs3L3fSIlI61157LePHj092GiISIdFFfw3Q3szamVltAoV9SeRKZtYQ6AW8FbasvpldXDAN9AM2JSRrERGRC0BCT++7+xkzmwC8D6QAc919s5mNC8bnBFcdBix392NhzVsAbwa7AdUEfufu7yUuexGJlfrpi1ROCR9wx92XAksjls2JmJ8PzI9Y9gXQuYLTE5E4UT99kcpHT+QTERGpJjS0rojE3bRp02jRokWy0xCRCCr6IhJ3d911V7JTEJEodHpfROLuiy++YO/e4nrjikiyxKXom1m/eLyOiFwY+vXrx9SpU5OdhohEKPPpfTN7DfiMwFP2/h5YHq+kREREJP7Kc03/OXf/CMDMlsUpHxG5AKifvkjlVKrT+2b2QMF0QcEPTv93PJMSkapP/fRFKp/SXtNPq5AsREREpMKV9vR+WzP7fviC4BP2RERCHnvsMfXTF6mESlv0FwLN+NsQubpwJyJFDB8+PNkpiEgUMRV9M/tHYCRwmMDIdv8NbHL3UxWXmohUVRs3bqR+/fqkpemKoEhlEus1/VeBd4BPCFzX/wmwuaKSEpGqbejQocyYMSPZaYhIhFhP7+9w9zeD04sqKhkRERGpOCUe6ZvZy2b2I+BjM/txYlISkapO/fRFKqfznd5/KbjOd4A7zex/zWyJmf3MzG6p+PREpKpSP32RyqfEou/uH7j7E+4+yt27AH8HPAL8D3BVWTZoZgPMbLuZ7TCzaVHivc3siJltCP78JNa2IiIiUrxSddlz9zME7twv0xP4zCwFmA1cD+QBa8xsibtviVj1I3cfXMa2IpJkTz/9NJdcckmy0xCRCOct+mZWB5gK9ACOAeuBhe6+swzbu4rATYFfBF/7dWAoEEvhLk9bEUmgIUOGJDsFEYkili57s4C2wEMEjrCbAvPM7I4ybK8VsCdsPi+4LFIPM/uzmS0zs/RStsXM7jWztWa2dv/+/WVIU0TKY/Xq1WzdujXZaYhIhFiKfmfgAXf/M3CVuz9A4BT7+DJsL9qdPZG3+a4H2rp7Z+AZ4A+laBtY6P68u2e7e3bz5s3LkKaIlMftt9/Or371q2SnISIRYj3SL1jvgJk9AYwCzpVhe3nApWHzrYG94Su4+1/d/WhweilQy8yaxdJWREREine+fvo9gHfd/ZvgoluAtUCD4HRprQHam1k7M6sNDAeWRGzzOxbs62NmVwVzPBhLWxGpHNRPX6RyOt+NfKOA2Wb2P8B7wHvu/ruybszdz5jZBOB9IAWY6+6bzWxcMD4HuBm4z8zOACeA4R74HyRq27LmIiIVS/30RSqfEou+u48DMLMrgIHAfDNrCHxI4EvAn9z9bGk2GDxlvzRi2Zyw6WeBZ2NtKyIiIrGJqZ++u28DtgFPmtlFwHUETu8/AWRXXHoiUhXNnTsX3UQrUvmUZ2jdBzS0rohE07dv32SnICJRaGhdEYm75cuXs2HDhmSnISIRNLSuiMTd3XffTb9+/fjtb3+b7FREJIyG1hUREakmNLSuiMSd+umLVE7n67L3AfBBwbyZ1QQ6Eng079XoVL+IFEP99EUqn/KMsvdKBecmIiIicRTLjXyzgHoERtkbzd9G2fuNuy+owNxEpIpatGgRTZo0SXYaIhIhlqLfGRjq7ofN7Cp37xk8+s8FVPRFpIgePXokOwURiSLRo+yJSDWwePFiPv3002SnISIRzlv03f2dOI6yJyLVwIQJE9RHX6QSivXhPAC4+2mgzKPsiYiISPLE+hheEZGYqZ++SOWkoi8iFUL99EUqn4QXfTMbYGbbzWyHmU2LEh9pZv8d/FltZp3DYrvNbKOZbTCztYnNXEREpGor1TX98jKzFGA2cD2QB6wxsyXuviVstV1AL3c/ZGYDgecJPP2vwHXufiBhSYtIqb3//vs0bNgw2WmISISEFn3gKgIj9n0BYGavA0OBUNF399Vh638CtE5ohiJSbpmZmclOQUSiSPTp/VbAnrD5vOCy4twNLAubd2C5ma0zs3srID8RiYOXX36ZVatWJTsNEYmQ6KIf7c6eqLf5mtl1BIr+g2GLe7p7V2AgMN7Mri2m7b1mttbM1u7fv7+8OYtIKU2dOpXf/U69e0Uqm0QX/Tzg0rD51sDeyJXMLBN4kcDjfw8WLHf3vcF/9wFvErhcUIS7P+/u2e6e3bx58zimLyIiUnUluuivAdqbWTszqw0MB5aEr2BmbYDFwJ3u/j9hy+ub2cUF00A/YFPCMheRmKmfvkjllNAb+dz9jJlNAN4HUoC57r7ZzMYF43OAnxAYye+5YD/fM+6eDbQA3gwuqwn8zt3fS2T+IhI79dMXqXwSffc+7r4UWBqxbE7Y9D3APVHafUFgxD8REREpg4QXfRG58H388cfUr18/2WmISAQVfRGJu7S0tGSnICJRqOiLSNzNnj2byy+/nL59+8a0fn5+Pnl5eZw8ebKCMxO5cNStW5fWrVtTq1atmNuo6ItI3M2YMYNbb7015qKfl5fHxRdfTGpqqm4AFImBu3Pw4EHy8vJo165dzO00yp6IJN3Jkydp2rSpCr5IjMyMpk2blvrsmIq+iMRdWfrpq+CLlE5Z/mZU9EWkQlyIRXz+/PlMmDABgDlz5vDyyy+XuP5PfvITVqxYUeI6ubm5rF69usR1KsKSJUuYOXMmAH/4wx/YsuVvg53Gkvf3v/99Dh8+zOHDh3nuuecqNNeyuOaaayp8G6mpqRw4UPygr7t37y70OOq1a9dy//33l2lbvXv3Zu3a8o8or2v6IiJlMG7cuPOu89Of/vS86+Tm5tKgQYOEFKlwQ4YMYciQIUCg6A8ePJiOHTsCseW9dGngcSu7d+/mueee41/+5V8qLtkySMYXqUgFRX/EiBEAZGdnk52dndScdKQvInG3adMmHnvssWSnUSo33ngjV155Jenp6Tz//POh5fPmzeOyyy6jV69e/OlPfwotnzFjBo8//jgAGzZsoHv37mRmZjJs2DAOHToEwOjRo3njjTeAwFHh9OnT6dq1KxkZGWzbto3du3czZ84cnnzySbp06cJHH31UKKeVK1fSpUsXunTpQlZWFt9++y0As2bNolu3bmRmZjJ9+nQgUGA6dOjA2LFjSU9Pp1+/fpw4cQKAnJwcOnbsSGZmJsOHDwf+dtZi9erVLFmyhClTptClSxd27twZynvZsmXceuutoXxyc3O54YYbQu/nwIEDTJs2jZ07d9KlSxemTJnCnXfeyVtvvRVqM3LkSJYsKfS09UKOHTvGoEGD6Ny5M506deLf//3fAVi3bh29evXiyiuvpH///nz11VdA4Ih30qRJXHvttXTo0IE1a9Zw00030b59ex555JHQ6zZo0CDq9t5++22uvvpqsrKy6Nu3L19//XXo93nXXXfRu3dv0tLSyMnJOe9no8Cjjz7K008/HZp/+OGHycnJYdq0aXz00Ud06dKFJ598ktzcXAYPHgzA0aNHGTNmDBkZGWRmZvL73/8egPvuu4/s7GzS09NDv9u4cvcL+ufKK690EanctmzZUmi+V69eRX5mz57t7u7Hjh2LGp83b567u+/fv79ILBYHDx50d/fjx497enq6HzhwwPfu3euXXnqp79u3z0+dOuXXXHONjx8/3t3dp0+f7rNmzXJ394yMDM/NzXV390cffdR/+MMfurv7qFGjfNGiRe7u3rZtW8/JyXF399mzZ/vdd99d5HUiDR482P/zP//T3d2//fZbz8/P9/fff9/Hjh3r586d87Nnz/qgQYN85cqVvmvXLk9JSfHPPvvM3d1vueUWf+WVV9zdvWXLln7y5El3dz906JC7u8+bNy/0XsLzDJ/Pz8/3Sy+91I8ePeru7uPGjQu9Ztu2bX3//v2+a9cuT09PD7XNzc31oUOHurv74cOHPTU11fPz84vd72+88Ybfc889ofnDhw/76dOnvUePHr5v3z53d3/99dd9zJgx7h74bEydOtXd3Z966ilv2bKl792710+ePOmtWrXyAwcOuLt7/fr1o27vm2++8XPnzrm7+wsvvOCTJ09298DvoUePHn7y5Enfv3+/N2nSxE+fPu3u0T8bkfsgKyvL3d3Pnj3raWlpfuDAAf/www990KBBoW2Hz0+dOjX0OSnIK3xbZ86c8V69evmf//zn0Ptes2ZNkfcT+bfj7g6s9WJqoo70RSTuZs6cybvvvpvsNEolJyeHzp070717d/bs2cPnn3/Op59+Su/evWnevDm1a9fmtttuK9LuyJEjHD58mF69egEwatQoVq1aFXUbN910EwBXXnklu3fvPm9OPXv2ZPLkyeTk5HD48GFq1qzJ8uXLWb58OVlZWXTt2pVt27bx+eefA9CuXTu6dOlSZBuZmZmMHDmSV199lZo1Y7+qW7NmTQYMGMDbb7/NmTNnePfddxk6dGiJbXr16sWOHTvYt28fr732Gj/4wQ9K3GZGRgYrVqzgwQcf5KOPPqJhw4Zs376dTZs2cf3119OlSxd+/vOfk5eXF2pTcFkiIyOD9PR0WrZsSZ06dUhLS2PPnj0l5peXl0f//v3JyMhg1qxZbN68ORQbNGgQderUoVmzZlxyySWhswDRPhvhUlNTadq0KZ999lnod9O0adMS81ixYgXjx48PzTdu3BiAhQsX0rVrV7Kysti8eXOhey3iQdf0RSTuZs2axYgRIxg0aFCZ2ufm5hYbq1evXonxZs2alRgvbnsrVqzg448/pl69evTu3TvUFSqeNyTWqVMHgJSUFM6cOXPe9adNm8agQYNYunQp3bt3Z8WKFbg7Dz30EP/8z/9caN3du3eHXr9gGwWn9999911WrVrFkiVL+NnPflao0J3PbbfdxuzZs2nSpAndunXj4osvPm+bO++8kwULFvD6668zd+7cEte97LLLWLduHUuXLuWhhx6iX79+DBs2jPT0dD7++OOobQreZ40aNQq95xo1ahTZrw8//HDoC+iGDRuYOHEikydPZsiQIeTm5jJjxowirwt/+x2V9NkId8899zB//nz+8pe/cNddd5W8gwicZY/8bO3atYvHH3+cNWvW0LhxY0aPHh33B1bpSF9Eqr0jR47QuHFj6tWrx7Zt2/jkk08AuPrqq8nNzeXgwYPk5+ezaNGiIm0bNmxI48aNQ9fjX3nlldBRfywuvvji0LX6SDt37iQjI4MHH3yQ7Oxstm3bRv/+/Zk7dy5Hjx4F4Msvv2Tfvn3Fvv65c+fYs2cP1113Hb/61a84fPhwqG0sOfTu3Zv169fzwgsvRD3TEa3t6NGjeeqppwBIT08P5dmnT58i7ffu3Uu9evW44447eOCBB1i/fj2XX345+/fvDxX9/Pz8Un1RCffYY4+xYcMGNmzYAAR+161atQLgpZdeOm/74j4bkYYNG8Z7773HmjVr6N+/P1Dyfu3Xrx/PPvtsaP7QoUP89a9/pX79+jRs2JCvv/6aZcuWleatxkRH+iISd16GfvrJNGDAAObMmUNmZiaXX3453bt3B6Bly5bMmDGDHj160LJlS7p27crZs2dD7QqO1F566SXGjRvH8ePHSUtLY968eTFv+4YbbuDmm2/mrbfe4plnnuEf/uEfQrGnnnqKDz/8kJSUFDp27MjAgQOpU6cOW7dupUePHkDghrVXX32VlJSUqK9/9uxZ7rjjDo4cOYK7M2nSJBo1alRoneHDhzN27FhycnJCNx4WSElJYfDgwcyfPz9qkWzatCk9e/akU6dODBw4kFmzZtGiRQs6dOjAjTfeGFrvq6++inqaf+PGjUyZMoUaNWpQq1Ytfv3rX1O7dm3eeOMN7r//fo4cOcKZM2f40Y9+FPoCUR4zZszglltuoVWrVnTv3p1du3aVuH5xn41ItWvX5rrrrqNRo0ah30VmZiY1a9akc+fOjB49mqysrND6jzzyCOPHj6dTp06kpKQwffp0brrpJrKyskhPTyctLY2ePXuW+/1Gsqr2x1la2dnZHo++jSISuyZNmnDHHXcUugO6JFu3bqVDhw4VnFV8TZw4ka5duzJmzJhkp1LpHD9+nIyMDNavX0/Dhg0BePbZZ2nTpk3oevyF5ty5c3Tt2pVFixbRvn37hG032t+Oma1z96h9A3V6X0SklB599FE+/fTTC7aAlceKFSu44oormDhxYqjgA0yYMOGC3V9btmzhe9/7Hn369ElowS+LhB/pm9kA4GkgBXjR3WdGxC0Y/z5wHBjt7utjaRuNjvRFEu/o0aPUqlWr0I1RJamKR/oilUGlPtI3sxRgNjAQ6AjcbmYdI1YbCLQP/twL/LoUbUWkEmjQoEHMBV9EEifRp/evAna4+xfufhp4HYjs9DkUeDn4jIFPgEZm1jLGtiJSCTz88MO8+eabpWpzod9fJBJvZfmbSXTRbwWEPzkhL7gslnViaVuhDh06ROvWrYv8PPPMMwDs2bMnarygn+rWrVujxhcuXAgEHjsZLV7wjOuVK1dGjRf0SV66dGnUeMHljYULF0aNb926FYC5c+dGjRc87OLZZ5+NGi945OjMmTOjxk+dOgUEroNGxsLHgZ40aVKReKdOnULxsWPHFomH30k7cuTIIvHw8dyHDh1aJB5+jbFv375F4gXPzAbo3r17kfjYsWND8U6dOhWJT5o0KRRPS0srEn/00UcBOH36dNR998tf/rJKfvZ+8YtfsHLlyqh/R9HUrVuXgwcPqvCLxMjdOXjwIHXr1i1Vu0R32Yv2lIvIv/Li1omlbeAFzO4lcGmANm3alCa/EtWqVYsBAwYUWZ6WlgbARRddFDXetm1bINBnM1q8oM9oo0aNosZbtGgBBB46Ei3erFmz0HrR4gXdc1q1ahU1XvCwjbZt20aNX3TRRUDgfUaL16pVCwg8ZCNavEaNwHfLDh06FImHdzPKyMgo0qe1fv36oemsrKwiRSH8qVfZ2dmhXAsU7FsIFO3mzZsXihf87iDw9LPU1NRC8YIBSCDwpLGDBw8Wiod3wenTpw/Hjh0rFM/IyAhN9+vXr8iDQwquxdWoUSPqvrvsssuAqvfZMzNGjRpVZP3itG7dmry8PPbv3x9zG5Hqrm7durRu3bpUbRJ6I5+Z9QBmuHv/4PxDAO7+y7B1fgPkuvtrwfntQG8g9Xxto9GNfCIiUp1Umhv5gDVAezNrZ2a1geFA5PBLS4B/soDuwBF3/yrGtiIiIlKMhJ7ed/czZjYBeJ9At7u57r7ZzMYF43OApQS66+0g0GVvTEltE5m/iIhIVaYn8omIiFxASjq9f8EXfTPbD/xvHF+yGXAgjq9XHWkflp/2YflpH8aH9mP5xXsftnX35tECF3zRjzczW1vcNyiJjfZh+Wkflp/2YXxoP5ZfIvehnr0vIiJSTajoi4iIVBMq+qX3fLITuABoH5af9mH5aR/Gh/Zj+SVsH+qavoiISDWhI30REZFqQkU/RmY2wMy2m9kOM5uW7HyqIjO71Mw+NLOtZrbZzH6Y7JyqKjNLMbPPzOydZOdSFZlZIzN7w8y2BT+PPZKdU1VjZpOCf8ebzOw1MyvdyC/VlJnNNbN9ZrYpbFkTM/sPM/s8+G/jitq+in4MzCwFmA0MBDoCt5tZx5JbSRRngB+7ewegOzBe+7HMfghsTXYSVdjTwHvufgXQGe3LUjGzVsD9QLa7dyLwlNThyc2qypgPRI5eNQ34wN3bAx8E5yuEin5srgJ2uPsX7n4aeB0YmuScqhx3/8rd1wenvyXwH21Ch0e+EJhZa2AQ8GKyc6mKzOz/AdcCvwVw99PufjipSVVNNYGLzKwmUA/Ym+R8qgR3XwV8E7F4KPBScPol4MaK2r6KfmxaAXvC5vNQsSoXM0sFsoBPk5xKVfQUMBU4l+Q8qqo0YD8wL3iJ5EUzq3++RvI37v4l8Djwf8BXBAZGW57crKq0FsGB5Qj+e0lFbUhFPzYWZZm6PZSRmTUAfg/8yN3/mux8qhIzGwzsc/d1yc6lCqsJdAV+7e5ZwDEq8HTqhSh4zXko0A74LlDfzO5IblYSCxX92OQBl4bNt0anssrEzGoRKPgL3H1xsvOpgnoCQ8xsN4HLTP9oZq8mN6UqJw/Ic/eCs0xvEPgSILHrC+xy9/3ung8sBq5Jck5V2ddm1hIg+O++itqQin5s1gDtzaydmdUmcMPKkiTnVOWYmRG4jrrV3Z9Idj5Vkbs/5O6t3T2VwOfwj+6uI6xScPe/AHvM7PLgoj7AliSmVBX9H9DdzOoF/677oJshy2MJMCo4PQp4q6I2VLOiXvhC4u5nzGwC8D6Bu1TnuvvmJKdVFfUE7gQ2mtmG4LJ/dfelyUtJqqmJwILgl/gvgDFJzqdKcfdPzewNYD2BXjmfoSfzxcTMXgN6A83MLA+YDswEFprZ3QS+UN1SYdvXE/lERESqB53eFxERqSZU9EVERKoJFX0REZFqQkVfRESkmlDRFxERqSZU9EVERKoJFX0REZFqQkVfRMrFzC4ys5VmlmJmqeHjhJfydWqb2argqG0iUgFU9EWkvO4CFrv72fK8SHDY6g+A2+KSlYgUoaIvIsUysw/N7Prg9M/NLCfKaiOJ8qxwM0sLDl3bLXgGYFtwGNtNZrbAzPqa2Z/M7HMzuyrY7A/B1xORCqCiLyIlmQ48bGYjgSxgUngw+Oz6NHffHbH8cgKjKY5x9zXBxd8DngYygSuAEcDfAw8A/xpcZxPQrULeiYhowB0RKZ67rwqOojYZ6B3lFH4z4HDEsuYEjvx/EDEw1S533whgZpuBD9zdzWwjkBrc3lkzO21mF7v7t/F/RyLVm470RaRYZpYBtAROFVOETwB1I5YdAfYQGFUx3Kmw6XNh8+cofABSBzhZ1pxFpHgq+iISlZm1BBYAQ4FjZtY/ch13PwSkmFl44T8N3Aj8k5mNKOU2mwL73T2/zImLSLFU9EWkCDOrBywGfuzuW4GfATOKWX05gWvzIe5+DBgMTDKzoaXY9HXA0lInLCIxMXdPdg4iUoWZWRYw2d3vjMNrLQYecvft5c9MRCLpSF9EysXdPwM+NLOU8rxOsCfAH1TwRSqOjvRFRESqCR3pi4iIVBMq+iIiItWEir6IiEg1oaIvIiJSTajoi4iIVBMq+iIiItWEir6IiEg18f8BpEsSmWiUk8wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dhdBC3h = np.array([M0.disvec(x,0) for x in X])[:,0].flatten()\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "f,s = plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., dhdBC3h, 'k--', mfc='none', label='adjoint sensitivity, semi-analytical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial h_{\\Gamma_1}$')\n",
    "plt.legend(loc=4);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.000000\n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f'% dhdBC3h[-1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c. Numerical"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAACdCAYAAACtmOGSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjS0lEQVR4nO3deXhU5d3/8fc3AYxBVBBEBEwIm+xJDBIIAgIqICJUFltksaBFZBG0GlCWioBPgapQ/aXUCvgAKiBWUdtSlZhCQQSEsgQkQFiUSkQTUUCW3L8/MuQhhGWAJCcz+byuyysz5z5zzmdORr65z5xz3+acQ0RERAJfiNcBREREpGCoqIuIiAQJFXUREZEgoaIuIiISJFTURUREgoSKuoiISJAo5XWAy1WxYkUXGRnpdQwREZEisXbt2m+dc5XO1hbwRT0yMpI1a9Z4HUNERKRImNnuc7Xp9LuIiEiQUFEXEREJEgF/+r2g3XPPPfz88895lnXu3Jlhw4aRnZ1Nx44d872me/fuPPTQQ/z000/cd999+dr79OlD7969ycjIoG/fvvnaH374Ybp168aePXsYNGhQvvZhw4bRoUMHvvzyS0aMGJGv/amnnqJVq1Zs2LCBp59+Ol/7uHHjaNq0KatWrWLixIn52idPnkzDhg1JTk7mD3/4Q772F154gZo1a/K3v/2NpKSkfO1JSUlUqVKFxYsXM2fOnHztc+bM4dprr2X+/PksWLAgX/ubb75JWFgYr732GkuWLMnTZmYsXrwYgJdffpmPP/44T3t4eDhz584FYNq0afz73//O037dddcxc+ZMACZNmsS6devytFetWpWXXnoJgLFjx5KampqnvVatWkyePBmAJ598kvT09DztDRs2ZOzYsUDO7+mbb77J0x4XF8dvf/tbAH7zm9+QlZWVp71ly5YMGTIEgH79+nHs2LE87e3bt2fAgAFkZ2fzwAMPcKbOnTvzq1/9isOHD/PQQw/la+/evTvdunXju+++Y9iwYfnae/fuTceOHfnqq69ITEzM1z5gwADatGnDzp07GT9+fJ626tWrM2TIEKpUqZLvdSLiDRX1M/z444/5ivrpz3/44Yd8rznV7pwjMzMzX/vRo0dz27/77rtztp88eZJvv/32nNs/fvw4Bw4cOGf7sWPH+O9//5uv/VShOHr0KF9//fU523/66Sf27duXr/348eNAzrHZs2fPOdt/+OGHs7afPHkSgMzMzHxFEXKOC8DBgwfZuXNnnjYzy32ckZFBWlpanvayZcvmPt6/fz9ffvllnvbrr78+9/HevXvZunVrnvbTf7fp6els2bIlT3tIyP+dzEpLS8v3+vDw8NzH27Zty/f+K1asmPt406ZNHDx4ME979erVcx+vX78+97NwSp06dXIfn+3akejoaCDnGK9evTpfe3x8PJDzPletWpWvvW3btgAcOXIk3x9EkPNHA8ChQ4dYvnx57nLnHHv37qVly5Yq6iLFiAX6hC5xcXFOF8qJFL09e/ZQvXp1zIwxY8awf/9+Ro0aRc2aNb2OJhLUzGytcy7ubG36Tl1ELslNN92UeyblxIkTzJ07l7p169K3b998ZzREpGioqIvIZZs8eTK7du1i+PDhvP3229SvX58XX3zR61giJY6KuogUiCpVqjBt2jTS09NJTEykTZs2AOzcuZO1a9d6G06khFBRF5ECValSJSZNmpR7Ed/kyZOJi4vj7rvvZuXKld6GEwlyKuoiUqimTZvGpEmTWL16NS1atKB9+/akpKR4HUskKKmoi0ihuvrqqxk1ahTp6elMnTqVTZs25Y494Jwj0O/AESlOVNRFpEiULVuWxx9/nF27duUOZLNs2TLi4+NZsmSJirtIAVBRF5EideWVV3LttdcCOQMiZWRk0KVLF2JiYli0aBHZ2dneBhQJYCrqIuKZTp06sW3bNmbPns2RI0fo0aMHd9xxh9exRAKWirqIeKp06dL069ePLVu28MYbb9C/f38gZ0Cb+fPn5w5DLCIXpqIuIsVCaGgo999/P3369AFgyZIl9O7dmzp16vCnP/0p35wMIpKfirqIFEtdu3bl/fffp3LlygwaNIiaNWsyffp09dxFzkNFXUSKJTPLHbBm6dKlREVFkZSURGhoKIAuqBM5CxV1ESnWzIw77riDlJQUli9fTkhICD/88AM333wzEydOzDdHvUhJpqIuIgGjQoUKAGRmZlKnTh2eeeYZIiIiGDt2LN99953H6US8p6IuIgHnpptu4v3332ft2rW0bduWCRMmEBERwVdffeV1NBFPeVbUzSzUzL4ws/d9zyuY2T/NbLvvZ3mvsolIYIiNjWXx4sVs3LiRxMREqlatCsDChQv5+uuvPU4nUvS87KkPB1JPe54IfOycqw187HsuInJBDRs25OmnnwYgKyuL/v37ExUVxaOPPsru3bs9TidSdDwp6mZWDbgbePW0xfcCc3yP5wBdiziWiASBa665ho0bN9K3b1/+/Oc/U6tWLQYOHMi+ffu8jiZS6Lzqqb8IPAmcfk9KZefcfgDfz+s9yCUiQSAqKoqZM2eyY8cOBg0axFtvvcWxY8eAnJHqRIJVkRd1M+sMHHDOrb2MbTxsZmvMbE1GRkYBphORYFK9enVmzJjB119/TVRUFAC/+MUv6NWrF//5z388TidS8LzoqScAXcwsHXgTaGtmc4FvzKwKgO/ngXNtwDk30zkX55yLq1SpUlFkFpEAVq5cOSBnwJpGjRrxt7/9jSZNmtC1a1fWrFnjcTqRglPkRd05N8o5V805FwncD3zinHsAeA/o51utH/BuUWcTkeAWEhLCxIkTSU9PZ/z48Xz66ac0bdqU119/3etoIgWiON2n/jxwh5ltB+7wPRcRKXAVKlRg3Lhx7N69m+eff57OnTsDkJKSQnJyMs45jxOKXBoL9A9vXFyc0+kzESkIHTp04B//+ActW7bkmWee4c4778TMvI4lkoeZrXXOxZ2trTj11EVEPPXOO+8wY8YM0tPT6dChA/Hx8XzyySdexxLx2yUVdTMra2ahBR1GRMRLV155JUOGDCEtLY2ZM2eSkZHBzp07ATh+/LhmhpNiz6+ibmYhZvYrM/vAzA4AW4H9ZrbZzKaYWe3CjSkiUnSuuOIKHnroIbZt20a/fjnX7yYlJdGoUSPmz5+ve92l2PK3p74MqAmMAm5wzlV3zl0P3AasAp43swcKKaOIiCdKly5N6dKlAahRowZmRu/evalXrx6zZs3i+PHjHicUycuvC+XMrLRz7ryfXn/WKQy6UE5Eikp2djbvvvsuEyZM4IsvvqBHjx4sWLDA61hSwlz2hXL+FGsvCrqISFEKCQmhW7durF27lg8++IDhw4cDcODAAaZPn87hw4c9Tigl3QWLupndYWZ/NrNo3/OHCz2ViEgxZmZ06tSJhIQEABYtWsTw4cOpUaMGU6ZM4ccff/Q4oZRU/vTUBwO/BR4ws7ZAdKEmEhEJMIMHD+bTTz+lSZMmPPnkk0RERDBp0iQNYiNFzp+inuGcy3TOPQHcCTQt5EwiIgGnVatWLF26lJUrV9KiRQvWr1+fO3DNTz/95HE6KSn8KeofnHrgnEsENEiyiMg5xMfHs2TJEubNmwdAamoqN954I0899RQHDpxzniqRAnHBou6ce/eM5zMKL46ISHA4dSvcFVdcQadOnZgyZQqRkZGMGDGCr7/+2uN0Eqwue5hYM1tREEFERIJRVFQUb7zxBqmpqfTs2ZMZM2bQsGFDXSkvhaIgxn6/sQC2ISIS1OrWrcvs2bP58ssvSUpKIjw8HOcckyZNIi0tzet4EiT8HSZ2hpk9bGbNzazcGc26vFNExE9RUVH07NkTgB07dvDss89St25d+vTpw9atWz1OJ4HO3576RqAxOXOcp5vZLjN7z8wmAmcWeRER8UOtWrXYtWsXI0aMYPHixdSvX59evXrx3//+1+toEqAuaT51M6tGTpFvBDRyznk27ruGiRWRYJCRkcELL7zAokWLWL9+PeHh4WRlZXHNNdd4HU2KmfMNE+vv2O/NgVWuGI6koKIuIsHk5MmThIaGcuLECW6++Wbq1KnDM888Q4sWLbyOJsXEZY/9DvQD1pnZm2bW38xuKLh4IiJySmhoKAAnTpxg4MCBfP755yQkJNCuXTuWLVumUerkvPyd0GWQcy4GGA+UB2ab2Uozm2RmrcwstDBDioiUNGFhYSQmJpKens60adPYsmULbdu25aOPPvI6mhRj/p5+L+WcO3HGsiuB24GOQPNznQoobDr9LiIlwdGjR3nrrbfo06cPISEhzJ49mwoVKnDPPffkDkcrJUNBnH5ffeYC59wR59yHzrmhXhV0EZGSIiwsjH79+hESEoJzjpdffpl7772XmJgYFi5cSHZ2ttcRpRjwt6jrz0ARkWLCzFi5ciVz5szhyJEj9OzZk4YNG7Js2TKvo4nHSvm5XiUzG3muRufcHwooj4iI+KFUqVL07duX3r17s3DhQiZOnEhYWBgAmZmZlC1bNnf8eSk5/O2phwJXkTPQzNn+ExERD4SGhnL//ffzn//8h+bNmwPw1FNPUbt2bZKSkvj55589TihFyd+e+n7n3LOFmkRERC7Z6RfLdevWjQ0bNvDII4/w3HPP8eSTTzJw4EDCw8M9TChFQd+pi4gEmQ4dOrBy5Ur++c9/UrNmTYYPH87o0aO9jiVFwN+eetsLrWBmVhxHnBMRKYnMjPbt29O+fXtSUlKIiIgAYPXq1SxdupShQ4dqCNog5G9P/R0zG2pmN52+0MzKmFlbM5tDzqhzIiJSzLRq1Sq3qC9dupQxY8YQERHB2LFjOXjwoMfppCD5W9Q7ACeBN8zsazPbYmY7ge3AL4EXnHOzCymjiIgUkGeeeYZ169bRrl07JkyYQGRkJJMmTfI6lhQQv06/O+eOAq8Ar5hZaaAicMQ5l1mI2UREpBDExMTw9ttvs3nzZiZOnJg7cE12djbffPMNVapU8TihXCp/e+q5nHPHnXP7L6egm1l1M1tmZqlmttnMhvuWVzCzf5rZdt/P8pe6DxEROb8GDRowf/58nn76aQDeeecdIiMjGTx4MLt37/Y4nVyKiy7qBeQE8Lhzrh4QDzxqZvWBROBj51xt4GPfcxERKUSnboe75ZZb6N+/P6+++iq1atViwIABpKWleZxOLoYnRd3X01/ne3wISAWqAvcCc3yrzQG6epFPRKQkioyM5E9/+hM7duxg0KBBzJs3j7vvvlvTvQYQr3rqucwsEogBPgMqO+f2Q07hB673MJqISIlUvXp1ZsyYwa5du5gzZw5mxpEjRxg4cCAbNmzwOp6cxwWLupmV9f28qqB37tvm28BjzrkfLuJ1D5vZGjNbk5GRUdCxREQEqFKlCvHx8QCsX7+eBQsWEB0dzb333oumvC6e/OmplzezIUDLgtyx7yr6t4F5zrnFvsXfmFkVX3sV4MDZXuucm+mci3POxVWqVKkgY4mIyFk0b96c3bt3M378eFJSUmjatCkdO3bk0KFDXkeT0/hT1NsB/YEoMzvv6XAze8qfnVrOVRl/AVLPmOHtPf5vEJt+wLv+bE9ERApf+fLlGTduHLt372by5MmUKVOGq67KOYm7c+dOffdeDNiFfglmVg8oDVR3zn1wRtuC058C0b4r1y+0zZbAv4CNQLZv8WhyvldfANwE7AF6OOe+O9+24uLinE4DiYh458CBA0RGRhIbG8uYMWO4884780wwIwXLzNY65+LO1nbBnrpzLtU5958zC7rPD865nr7/egAf+RPIObfcOWfOucbOuWjffx865w4659o552r7fp63oIuIiPeuvvpqpkyZwp49e+jQoQPNmjXjvffeU8/dAxd19buZvWFmT5lZRzOrBDx3xipPF1w0EREJBGFhYTz66KOkpaXx5z//mYMHD9K1a1e2b9/udbQS52JvaUsCDgP3Af8EXjCzcqca1bMWESm5ypQpw8CBA9m2bRuffPIJderUASAxMZG5c+dy4sQJjxMGv4st6jcBnwKDnHPRwCLgdwUdSkREAlepUqVo06YNAEePHuXvf/87ffr0oV69erz22mscP37c24BB7GKLeh1yTrlvNbMvgDuBX5hZO9/peBERkVxhYWGsW7eOd955h6uvvpoBAwZQu3ZtVq5c6XW0oHSxRX26c66Lc64W0BaYC4QDvYEPCzqciIgEvpCQELp27cqaNWv44IMPiIqKIioqCoBdu3Zx+PBhjxMGjwve0pZn5ZzeeXlyxmrfBtwMbHPODS+ceBemW9pERAKTc46EhAR27NjBE088wSOPPJJ737uc22Xd0nY651wMUBN4AkgB/gg8ftkJRUSkxDEznn/+eZo0acKTTz5JREQEzz33HFlZWV5HC1iXMp/6SefcZufcYufc+845Xc4oIiKXpFWrVixdupRVq1bRokULxowZw+uvv+51rIBVyusAIiIizZo1Y8mSJXzxxRfUrVsXgLlz57Jx40Yef/xxrr9ek3b6w/OpV0VERE6JiYkhPDwcgI0bNzJ16lQiIyMZMWIEX331lcfpij8VdRERKZb+53/+h9TUVHr16sWMGTOIiopi6tSpXscq1lTURUSk2KpTpw6zZs1i+/btPPjgg9SqVQuAzMxM0tLSPE5X/Kioi4hIsVejRg2SkpLo2rUrANOnT6du3bo88MADpKamehuuGFFRFxGRgPPwww8zcuRI/vrXv9KgQQN69uzJhg0bvI7lORV1EREJODfccANTpkwhPT2dUaNG8fe//53Ro0d7HctzFzWiXHF0thHljh8/zr59+zh69KhHqUTOLSwsjGrVqlG6dGmvo4gEje+//56srCwiIyPZuXMnQ4cOZfTo0SQkJHgdrcCdb0S5oLxPfd++fZQrV47IyEjMzOs4Irmccxw8eJB9+/ZRo0YNr+OIBI3y5ctTvnx5AL788ktWr15Ny5Ytuf322xkzZgxt2rQpEfUgKE+/Hz16lOuuu65E/AIlsJgZ1113nc4iiRSiDh06kJ6ezrRp00hNTaVt27a0a9eOkydPeh2t0AVlUQdU0KXY0mdTpPCVLVuWkSNHsmvXLv74xz+SkJBAaGgoAMuXLyfQv3o+l6At6oFg9uzZDBkyBICkpKQLjnc8duxYPvroo/Ouk5yczL///e/Lzpaenk7Dhg0vuN6kSZPyPG/RosVl77s46tSpE5mZmRf9utN/xyJS9MLCwnj00UeZMGECAKtXr+a2224jJiaGhQsXkp2d7XHCgqWiXkwMGjSIvn37nnedZ599lvbt2593nYIq6v46s6gX5b6LgnOO7OxsPvzwQ6699lqv44jIZYqNjWXOnDkcPXqUnj170rBhQ+bNm8eJE8ExN5mKeiHp2rUrt9xyCw0aNGDmzJm5y2fNmkWdOnVo3bo1K1asyF0+fvz43OEP169fT3x8PI0bN6Zbt258//33APTv359FixYBEBkZybhx44iNjaVRo0Zs3bqV9PR0kpKSeOGFF4iOjuZf//pXnkyffvop0dHRREdHExMTw6FDh3DO8dvf/paGDRvSqFEj3nrrrXzv5czeZufOnUlOTiYxMZEjR44QHR1N7969AXLnQj7XdpOTk2nTpg3du3fn5ptvpnfv3hc8Dda/f3+GDRtGixYtiIqKyj0GycnJdO7cOXe9IUOGMHv27NzjM3r0aJo3b05cXBzr1q3jrrvuombNmiQlJeW+ZsqUKTRt2pTGjRszbtw4IOcsRb169Rg8eDCxsbHs3buXyMhIvv32WwBef/11GjduTJMmTejTpw8AS5YsoVmzZsTExNC+fXu++eab874nEfFGqVKl6Nu3L5s3b+bNN98kNDSUoUOHcvjwYa+jFYigvPr9TG3atMm3rGfPngwePJjDhw/TqVOnfO39+/enf//+fPvtt3Tv3j1PW3Jy8gX3+dprr1GhQgWOHDlC06ZNue+++zh27Bjjxo1j7dq1XHPNNdx+++3ExMTke23fvn2ZMWMGrVu3ZuzYsfzud7/jxRdfzLdexYoVWbduHa+88gpTp07l1VdfZdCgQVx11VU88cQT+dafOnUqL7/8MgkJCfz444+EhYWxePFi1q9fz4YNG/j2229p2rQprVq1uuD7A3j++ef54x//yPr16/O1nW+7X3zxBZs3b+bGG28kISGBFStW0LJly/Pua//+/SxfvpytW7fSpUuXfL+Ts6levTorV65kxIgR9O/fnxUrVnD06FEaNGjAoEGDWLp0Kdu3b2f16tU45+jSpQspKSncdNNNbNu2jVmzZvHKK6/k2ebmzZuZOHEiK1asoGLFinz33XcAtGzZklWrVmFmvPrqq/z+979n2rRpfh1HESl6oaGh9OrVix49epCWlsbVV19NdnY299xzD507d+bXv/41V1xxhdcxL5p66oVk+vTpNGnShPj4ePbu3cv27dv57LPPaNOmDZUqVaJMmTL06tUr3+uysrLIzMykdevWAPTr14+UlJSz7uMXv/gFALfccgvp6ekXzJSQkMDIkSOZPn06mZmZlCpViuXLl/PLX/6S0NBQKleuTOvWrfn8888v/Y37nG+7t956K9WqVSMkJITo6Gi/snft2pWQkBDq16/vdy+4S5cuADRq1IhmzZpRrlw5KlWqRFhYGJmZmSxdupSlS5cSExNDbGwsW7duZfv27QBEREQQHx+fb5uffPIJ3bt3p2LFigBUqFAByLmN8q677qJRo0ZMmTKFzZs3+5VRRLwVEhJCnTp1AMjIyOD7779n8ODBREVF8dJLLwVcD75E9NTP17MODw8/b3vFihX96pmfub+PPvqIlStXEh4eTps2bXJvYSrIK59P/RUZGhrq1/dBiYmJ3H333Xz44YfEx8fz0Ucf+XUFaKlSpfJcTOLP7Vjn2+7pf/36m/3015za9oVynXpNSEhInteHhIRw4sQJnHOMGjWK3/zmN3lel56eTtmyZc+awzl31t/h0KFDGTlyJF26dCE5OZnx48df8D2JSPFSuXJlVqxYwSeffMKECRN47LHHmDRpEsuWLaN+/fpex/OLeuqFICsri/LlyxMeHs7WrVtZtWoVAM2aNSM5OZmDBw9y/PhxFi5cmO+111xzDeXLl8/9Pvx///d/c3vt/ihXrhyHDh06a9uOHTto1KgRTz31FHFxcWzdupVWrVrx1ltvcfLkSTIyMkhJSeHWW2/N87rIyEjWr19PdnY2e/fuZfXq1bltpUuX5vjx4/n25c92zzRq1Cjeeecdv99rREQEW7Zs4eeffyYrK4uPP/7Y79cC3HXXXbz22mv8+OOPAHz11VccOHDgvK9p164dCxYs4ODBgwC5p9+zsrKoWrUqAHPmzLmoHCJSfJgZ7dq1Izk5mZSUFDp16pTbk1+5cuUl3QVTlEpET72odejQgaSkJBo3bkzdunVzT+NWqVKF8ePH07x5c6pUqUJsbGyewRBO9QDnzJnDoEGDOHz4MFFRUcyaNcvvfd9zzz10796dd999lxkzZnDbbbfltr344ossW7aM0NBQ6tevT8eOHSlTpgwrV66kSZMmmBm///3vueGGG/KcEk9ISKBGjRo0atSIhg0bEhsbm9v28MMP07hxY2JjY5k3b17u8m7dup11u1u3bj1n9o0bN+aeMvdH9erV6dmzJ40bN6Z27dpnvT7hfO68805SU1Np3rw5kHOR39y5c3PvZT2bBg0a8PTTT9O6dWtCQ0OJiYlh9uzZjB8/nh49elC1alXi4+PZtWvXRWURkeLntttuy/039NixY3Tr1o0jR44wbNgwHnvsMa677jqPE+YXlGO/p6amUq9ePY8SXZqhQ4cSGxvLgw8+6HUUz9x111384x//8DpGkQjEz6hISbd+/Xqee+453n77bcqWLcvgwYN5/PHHqVy5cpHmON/Y7zr9XgyMGTOGzz777KJ6qcGopBR0EQlM0dHRLFq0iE2bNtGlSxemTZvGpk2bvI6Vh3rqIh7QZ1Qk8KWnpxMREYGZkZiYSFZWFomJiURERBTqftVTFxERKWCnzwT6888/85e//IVatWoxYMAA0tLSPMlU7Iq6mXUws21mlmZmiZe6nUA/AyHBS59NkeDzwgsvsHPnTh555BHmz59P3bp184xeWVSKVVE3s1DgZaAjUB/4pZld9M2BYWFhHDx4UP94SrFzaj71sLAwr6OISAGrVq0a06dPZ9euXYwcOdLv0TkLUnG7pe1WIM05txPAzN4E7gW2XMxGqlWrxr59+8jIyCiEiCKXJywsjGrVqnkdQ0QKyQ033MCUKVM82XdxK+pVgb2nPd8HNLvYjZQuXZoaNWoUWCgREZFAUKxOvwNnG0M13zl0M3vYzNaY2Rr1xkVERHIUt6K+D6h+2vNqwNdnruScm+mci3POxVWqVKnIwomIiBRnxa2ofw7UNrMaZlYGuB94z+NMIiIiAaHYDT5jZp2AF4FQ4DXn3MQLrJ8B7C7ACBWBbwtweyWRjuHl0zEsGDqOl0/H8PIV9DGMcM6d9TR1sSvqXjOzNecaqUf8o2N4+XQMC4aO4+XTMbx8RXkMi9vpdxEREblEKuoiIiJBQkU9v5leBwgCOoaXT8ewYOg4Xj4dw8tXZMdQ36mLiIgECfXURUREgoSKuk9BzQ5XkplZdTNbZmapZrbZzIZ7nSlQmVmomX1hZu97nSUQmdm1ZrbIzLb6Po/Nvc4UaMxshO//401m9oaZaRYiP5jZa2Z2wMw2nbasgpn908y2+36WL6z9q6hTcLPDCSeAx51z9YB44FEdx0s2HEj1OkQAewn4u3PuZqAJOpYXxcyqAsOAOOdcQ3LGDbnf21QBYzbQ4YxlicDHzrnawMe+54VCRT1H7uxwzrljwKnZ4eQiOOf2O+fW+R4fIucf0qrepgo8ZlYNuBt41essgcjMrgZaAX8BcM4dc85lehoqMJUCrjSzUkA4ZxmyW/JzzqUA352x+F5gju/xHKBrYe1fRT3H2WaHUzG6DGYWCcQAn3kcJRC9CDwJZHucI1BFARnALN9XGK+aWVmvQwUS59xXwFRgD7AfyHLOLfU2VUCr7JzbDzmdH+D6wtqRinoOv2aHE/+Y2VXA28BjzrkfvM4TSMysM3DAObfW6ywBrBQQC/w/51wM8BOFeLozGPm+870XqAHcCJQ1swe8TSX+UFHP4dfscHJhZlaanII+zzm32Os8ASgB6GJm6eR8DdTWzOZ6Gyng7AP2OedOnSVaRE6RF/+1B3Y55zKcc8eBxUALjzMFsm/MrAqA7+eBwtqRinoOzQ5XAMzMyPkeM9U59wev8wQi59wo51w151wkOZ/DT5xz6iFdBOfcf4G9ZlbXt6gdsMXDSIFoDxBvZuG+/6/boYsNL8d7QD/f437Au4W1o1KFteFA4pw7YWZDgH/wf7PDbfY4ViBKAPoAG81svW/ZaOfch95FkhJqKDDP90f6TuBBj/MEFOfcZ2a2CFhHzl0tX6CR5fxiZm8AbYCKZrYPGAc8DywwswHk/MHUo9D2rxHlREREgoNOv4uIiAQJFXUREZEgoaIuIiISJFTURUREgoSKuoiISJBQURcREQkSKuoiIiJBQkVdRM7LzK40s099c7xHnj5P9EVup4yZpfhm/RKRQqCiLiIX8mtgsXPu5OVsxDet8cdArwJJJSL5qKiLlGBmtszM7vA9fs7Mpp9ltd6cZaxqM4vyTW3a1NeD3+qb5nSTmc0zs/ZmtsLMtpvZrb6X/dW3PREpBCrqIiXbOOBpM+sNxAAjTm/0jZ0e5ZxLP2N5XXJm43vQOfe5b3Et4CWgMXAz8CugJfAEMNq3ziagaaG8ExHRhC4iJZlzLsU3C9dIoM1ZTrFXBDLPWFaJnJ77fWdMfLTLObcRwMw2Ax8755yZbQQiffs7aWbHzKycc+5Qwb8jkZJNPXWREszMGgFVgJ/PUWSPAGFnLMsC9pIzK9/pfj7tcfZpz7PJ24G4Ajh6qZlF5NxU1EVKKDOrAswD7gV+MrO7zlzHOfc9EGpmpxf2Y0BXoK+Z/eoi93kdkOGcO37JwUXknFTURUogMwsHFgOPO+dSgQnA+HOsvpSc78ZzOed+AjoDI8zs3ovY9e3AhxcdWET8ovnUReS8zCwGGOmc61MA21oMjHLObbv8ZCJyJvXUReS8nHNfAMvMLPRytuO7kv6vKugihUc9dRERkSChnrqIiEiQUFEXEREJEirqIiIiQUJFXUREJEioqIuIiAQJFXUREZEgoaIuIiISJP4/asBLxwhN4L8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "M0 = flopy.modflow.Modflow(modelname='model', exe_name='../mf2005.exe')\n",
    "flopy.modflow.ModflowDis(M0, nlay=1, nrow=1, ncol=ncol, nper=1, delr=1., delc=1., top=0., botm=-b, steady=True, \n",
    "                         perlen=1., nstp=1)\n",
    "flopy.modflow.ModflowBas(M0, ibound=np.ones([nrow, ncol]), strt=BC3h*np.ones([nrow, ncol], dtype=float))\n",
    "flopy.modflow.ModflowLpf(M0, hk=K, vka=-999., ss=-999., sy=-999., ipakcb=53)\n",
    "flopy.modflow.ModflowGhb(M0, stress_period_data={0: [0, 0, ncol-1, 0., BC3c]})\n",
    "flopy.modflow.ModflowWel(M0, stress_period_data={0: [0, 0, ocol, 1.]})\n",
    "flopy.modflow.ModflowPcg(M0, hclose=1e-6, rclose=1e-6)\n",
    "flopy.modflow.ModflowOc(M0, stress_period_data={(0,0): ['save head', 'save budget']})\n",
    "M0.write_input()\n",
    "success, buff = M0.run_model(silent=True)\n",
    "A = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X/1000., A, 'k--', mfc='none', label='adjoint solution, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\psi^*_1$ (T/$L^2$)')\n",
    "plt.legend(loc=3);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgsAAACdCAYAAAA6/zu5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcg0lEQVR4nO3deZgU1b3/8fdHMCIRN9xQlhlzXQAHBhkJYAygRBFRQdkSJJB4XSJ6xe2KelGfmF/iDRqNPxWDcYuSGNxRMSrEQeMWQYkRARUlYZTIoqiIIgPf+0cXk56tp5HpmWb4vJ6nn64651T1t4t66O+cOlVHEYGZmZlZbbZr7ADMzMwsvzlZMDMzs4ycLJiZmVlGThbMzMwsIycLZmZmlpGTBTMzM8uoeWMHkK/22GOPKCgoaOwwzMzMGszcuXNXRsSeVcvzJlmQtDvwR6AAWAIMj4iPa2h3LnAaIODWiLg+Kb8yKV+RNL00ImYkdZcApwIbgP+KiCfriqegoIA5c+Zs0XcyMzPbmkj6R03l+XQZYgIwKyIOAGYl65VIOoRUQtAD6AoMknRAWpPrIqI4eW1KFDoBI4HOwADgZknNcvtVzMzMmo58ShZOBO5Klu8CBtfQpiPwUkSsjYhyYDYwJIv93hsR6yLiPeAdUsmGmZmZZSGfkoW9I2IZQPK+Vw1t3gC+K6m1pJbAQKBdWv3Zkl6XdLuk3ZKy/YClaW3KkjIzMzPLQoOOWZA0E9inhqrLstk+IhZI+l/gaWAN8DegPKmeDFwFRPJ+LfBjUmMbqu2qlvhOB04HaN++fTYhmdk2Yv369ZSVlfHll182dihmW6xFixa0bduW7bffPqv2DZosRET/2uokfSipTUQsk9QGWF7LPm4Dbku2+TmpngIi4sO0fd0KPJasllG596Et8EEt+54CTAEoKSnxDFtmVqGsrIxWrVpRUFCAVNPfIGZbh4hg1apVlJWVUVhYmNU2+XQZYjowJlkeAzxSUyNJeyXv7YGTgD8k623Smg0hdcli035HStpBUiFwAPDXeo/ezJq0L7/8ktatWztRsK2eJFq3br1ZvWR5c+skcDUwTdKpwD+BYQCS9gV+GxEDk3YPSGoNrAfGpd1e+UtJxaQuMSwBzgCIiPmSpgFvkrpkMS4iNjTMVzKzpsSJgjUVm3su503PQkSsioijIuKA5P2jpPyDtESBiDgiIjpFRNeImJVWPjoiiiKiS0ScsGmwZFL3/yLiWxFxUEQ80bDfzMys4dx5552cffbZANxyyy387ne/y9j+8ssvZ+bMmRnblJaW8sILL9RbjNmaPn06V199NQAPP/wwb775ZkVdNnEPHDiQ1atXs3r1am6++eacxppr6cdicxUUFLBy5cot+vx86lkwM7N6dOaZZ9bZ5qc//WmdbUpLS9lpp53o3bt3fYSVtRNOOIETTjgBSCULgwYNolOnTkB2cc+YMQOAJUuWcPPNN3PWWWflLtgcKi8vr3QsGkPe9CyYmVlmgwcPpnv37nTu3JkpU6ZUlN9xxx0ceOCB9OnTh+eff76i/Morr+Saa64BYN68efTs2ZMuXbowZMgQPv44dQV37Nix3H///UDqL9ArrriCQw89lKKiIhYuXMiSJUu45ZZbuO666yguLua5556rFNPs2bMpLi6muLiYbt268dlnnwEwadIkDjvsMLp06cIVV1wBpH60O3bsyGmnnUbnzp05+uij+eKLLwC44YYb6NSpE126dGHkyJHAv3tJXnjhBaZPn85FF11EcXExixcvroj7iSeeYPjw4RXxlJaWcvzxx1d8n5UrVzJhwgQWL15McXExF110EaNHj+aRR/49LG7UqFFMnz691uOeKe6+fftWPO135cqVbJom4M4772Tw4MEcf/zxFBYWcuONN/KrX/2Kbt260bNnTz766CMAFi9ezIABA+jevTtHHHEECxcurPh3Of/88+nXrx8XX3xxpR6jDz/8kCFDhtC1a1e6du1a0etT2/lRLyLCrxpe3bt3DzOzTd58883GDiFWrVoVERFr166Nzp07x8qVK+ODDz6Idu3axfLly2PdunXRu3fvGDduXEREXHHFFTFp0qSIiCgqKorS0tKIiJg4cWKce+65ERExZsyYuO+++yIiokOHDnHDDTdERMRNN90Up556arX9VDVo0KD4y1/+EhERn332Waxfvz6efPLJOO2002Ljxo2xYcOGOO6442L27Nnx3nvvRbNmzeK1116LiIhhw4bF3XffHRERbdq0iS+//DIiIj7++OOIiLjjjjsqvkt6nOnr69evj3bt2sWaNWsiIuLMM8+s2GeHDh1ixYoV8d5770Xnzp0rti0tLY0TTzwxIiJWr14dBQUFsX79+lqPe6a4+/TpE6+88kpERKxYsSI6dOhQEfu3vvWt+PTTT2P58uWx8847x+TJkyMiYvz48XHddddFRMSRRx4Zb731VkREvPTSS9GvX7+K73fcccdFeXl5tWMxfPjwiu3Ly8tj9erVEVHz+ZF+HKqq6ZwG5kQNv4m+DGFm9jX07du3Wtnw4cM566yzWLt2LQMHDqxWP3bsWMaOHcvKlSsZOnRopbrS0tI6P/OGG27goYceAmDp0qW8/fbb/Otf/6Jv377suWdq7p8RI0bw1ltvVdruk08+YfXq1fTp0weAMWPGMGzYsBo/46STTgKge/fuPPjgg3XGdPjhh3P++eczatQoTjrpJNq2bctTTz3FU089Rbdu3QBYs2YNb7/9Nu3bt6ewsJDi4uKKz1iyZAkAXbp0YdSoUQwePJjBgwfX+bmbNG/enAEDBvDoo48ydOhQHn/8cX75y19m3KZPnz6MGzeO5cuX8+CDD3LyySfTvHnmn8Pa4s6kX79+tGrVilatWrHLLrtU9HgUFRXx+uuvs2bNGl544YVK/xbr1q2rWB42bBjNmlWfneDPf/5zxViUZs2ascsuuwA1nx+tW7euM85sOFkwM9sKlJaWMnPmTF588UVatmxJ3759K259q8+7NHbYYQcg9SNUXl5eR2uYMGECxx13HDNmzKBnz57MnDmTiOCSSy7hjDPOqNR2yZIlFfvf9BmbuvMff/xxnn32WaZPn85VV13F/Pnzs455xIgR3HTTTey+++4cdthhtGrVqs5tRo8ezdSpU7n33nu5/fbb62xfW9zNmzdn48aNANVuRUzfZrvttqtY32677SgvL2fjxo3suuuuzJs3r8bP/OY3v1lnXJtkOj/qg5MFM7OvIVNPQMuWLTPW77HHHln1JKT75JNP2G233WjZsiULFy7kpZdeAuDb3/425557LqtWrWLnnXfmvvvuo2vXrpW23WWXXdhtt9147rnnOOKII7j77rsrehmy0apVKz799NMa6xYvXkxRURFFRUW8+OKLLFy4kGOOOYaJEycyatQodtppJ95///2MTwrcuHEjS5cupV+/fnznO9/h97//PWvWrKkWw6bxEFX17duXU089lVtvvZURI0bUGH/VbceOHUuPHj3YZ5996Ny5MwDvv/8+P/zhD5k1a1a1fdSmoKCAuXPn0qNHj4qxH9naeeedKSws5L777mPYsGFEBK+//nq1f7+qjjrqKCZPnsz48ePZsGEDn3/+ea3nR33xAEczs63AgAEDKC8vp0uXLkycOJGePXsC0KZNG6688kp69epF//79OfTQQyttt6nX4a677uKiiy6iS5cuzJs3j8svvzzrzz7++ON56KGHahzgeP3113PIIYfQtWtXdtxxR4499liOPvpofvCDH9CrVy+KiooYOnRorT/0ABs2bOCUU06hqKiIbt26cd5557HrrrtWajNy5EgmTZpEt27dWLx4caW6Zs2aMWjQIJ544gkGDRpUbf+tW7fm8MMP55BDDuGiiy4CYO+996Zjx4786Ec/qmi3bNmyOi9HVHXhhRcyefJkevfu/bVuT5w6dSq33XYbXbt2pXPnzpUGXtbm17/+Nc888wxFRUV0796d+fPn13p+1BelxjNYVSUlJbFphKuZ2YIFC+jYsWNjh7FZzjnnHA499NBKP4iWsnbtWoqKinj11VcrrvnfeOONtG/fvlFvUWxINZ3TkuZGREnVtu5ZMDNrgiZOnMjLL7+8zfzwbY6ZM2dy8MEHc84551QkCgBnn322j1ct3LNQC/csmFm6rbFnwSwT9yyYmZlZvXGyYGaWJffEWlOxueeykwUzsyy0aNGCVatWOWGwrV5EsGrVKlq0aJH1Nn7OgplZFtq2bUtZWRkrVqxo7FDMtliLFi1o27Zt1u2dLJiZZWH77bensLCwscMwaxS+DGFmZmYZ5U2yIGl3SU9Lejt5362WdudKekPSfEnj08r/KGle8loiaV5SXiDpi7S6WxrmG5mZmTUN+XQZYgIwKyKuljQhWb84vYGkQ4DTgB7AV8CfJD0eEW9HxIi0dtcCn6RtujgiinP9BczMzJqivOlZAE4E7kqW7wIG19CmI/BSRKyNiHJgNjAkvYFSD0IfDvwhd6GamZltO/IpWdg7IpYBJO971dDmDeC7klpLagkMBNpVaXME8GFEvJ1WVijpNUmzJR2Ri+DNzMyaqga9DCFpJrBPDVWXZbN9RCyQ9L/A08Aa4G9A1QnXv0/lXoVlQPuIWCWpO/CwpM4RUW2+VUmnA6cDtG/fPpuQzMzMmrwGTRYion9tdZI+lNQmIpZJagMsr2UftwG3Jdv8HChL20dz4CSge1r7dcC6ZHmupMXAgUC1iR8iYgowBVJzQ2z2FzQzM2uC8ukyxHRgTLI8BqhxUm9JeyXv7UklBum9CP2BhRGRnkDsKalZsrw/cADwbr1Hb2Zm1kTl090QVwPTJJ0K/BMYBiBpX+C3ETEwafeApNbAemBcRHycto+RVB/Y+F3gp5LKgQ3AmRHxUQ6/h5mZWZPiKapr4SmqzcxsW+Mpqs3MzOxrcbJgZmZmGTlZMDMzs4y+9gBHST9MFncG1kXErfUTkpmZmeWTLelZ2Bs4OCJuBPatp3jMzMwsz2zJrZOzI+KvyfKr9RGMmZmZ5Z/N6lmQdOGm5bREgYh4tD6DMjMzs/yxuZch9s9JFGZmZpa3NvcyRAdJA9MLImJGPcZjZmZmeWZzk4VpwJ5p6378o5mZWROXMVmQNBFYGxHXAkTEXQ0SlZmZmeWNunoWRgPFVQsl/SewZ0T8IhdBmZmZWf6oa4DjFxGxtobyu4FTchCPmZmZ5Zk6kwVJbaoWRsQ6oDw3IZmZmVk+qStZuBZ4RFKH9EJJewEbcxaVmZmZ5Y2MYxYi4j5JLYG5kl4C5pFKMIYBV+Y8OjMzM2t0dT6UKbkDopDUbZPbA18C34+IqfUZiKRhkuZL2iipJEO7AZIWSXpH0oS08t0lPS3p7eR9t7S6S5L2iyQdU59xm5mZNXVZPWchIj4DfpfjWN4ATgJ+U1sDSc2Am4DvAWXAK5KmR8SbwARgVkRcnSQRE4CLJXUCRgKdSU14NVPSgRGxIbdfx8zMrGmoM1mQtAPw30Av4HNSk0ZNi4jF9RlIRCxIPi9Tsx7AOxHxbtL2XuBE4M3kvW/S7i6gFLg4Kb83GZT5nqR3kv28WJ/xZ6Nv377VyoYPH85ZZ53F2rVrGThwYLX6sWPHMnbsWFauXMnQoUOr1f/kJz9hxIgRLF26lNGjR1erv+CCCzj++ONZtGgRZ5xxRrX6//mf/6F///7MmzeP8ePHV6v/+c9/Tu/evXnhhRe49NJLq9Vff/31FBcXM3PmTH72s59Vq//Nb37DQQcdxKOPPsq1115brf7uu++mXbt2/PGPf2Ty5MnV6u+//3722GMP7rzzTu68885q9TNmzKBly5bcfPPNTJs2rVp9aWkpANdccw2PPfZYpbodd9yRJ554AoCrrrqKWbNmVapv3bo1DzzwAACXXHIJL75Y+ZRp27Yt99xzDwDjx49n3rx5leoPPPBApkyZAsDpp5/OW2+9Vam+uLiY66+/HoBTTjmFsrKySvW9evXiF79I3Z188skns2rVqkr1Rx11FBMnTgTg2GOP5YsvvqhUP2jQIC68MDWdi889n3vpfO41jXOvIWUzN8QkoANwCam/5lsDd0hqjFsn9wOWpq2XJWUAe0fEMoDkfa8stqlE0umS5kias2LFinoN3MzMbGuliMxPbJY0GzgxIlZLej4iDk96G0ojotdmfZg0E9inhqrLIuKRpE0pcGFEzKlh+2HAMRHxn8n6aKBHRJwjaXVE7JrW9uOI2E3STcCLEXFPUn4bMCMiHsgUa0lJScyZUy0EMzOzJkvS3IioNm4wmzELk/h3D8RKSb8CFvI1bp2MiP6bu00VZUC7tPW2wAfJ8oeS2kTEsuTZEMuz2MbMzMzqkPEyhKRewOMR8VFSNAyYA+yULDe0V4ADJBVK+gapgYvTk7rpwJhkeQzwSFr5SEk7SCoEDgD+2oAxm5mZbdXq6lkYA9wk6S3gT8CfIuL3uQhE0hDg/5Oa1fJxSfMi4hhJ+wK/jYiBEVEu6WzgSaAZcHtEzE92cTUwTdKpwD9JkpmImC9pGqlBkOXAON8JYWZmlr06xywASDoYOBY4BtgFeIZU8vB8U/3h9ZgFMzPb1tQ2ZiGbuyGIiIURcV1EDACOBP5C6i/3l+s3TDMzM8s3WT2USdKRwChgNamHJ71O6o6FdbkLzczMzPJBVskCcA8wLmnfBRhM6omI/5GbsMzMzCxfZJssvBMRDyXL9+UqGDMzM8s/dd06+TtJ44EXJV3QMCGZmZlZPqlrgONdSZt9gNGS/iFpuqSrkqcpmpmZWROX8TJERMwCKmY4kdQc6AR0Bb6NL0mYmZk1eVsy6+TdOY7NzMzM8sCWzDo5KpeBmZmZWX7I5m6Irvx71ske6bNOAlNzGp2ZmZk1ugadddLMzMy2PnVehoiIx/Jo1kkzMzNrYNk+lAmAiPgKyMmsk2ZmZpafsppIyszMzLZdThbMzMwso7xJFiQNkzRf0kZJ1ebSTms3QNIiSe9ImpBWPknSQkmvS3pI0q5JeYGkLyTNS163NMDXMTMzazLyJlkgNfX1ScCztTWQ1Ay4CTiW1JMkvy+pU1L9NHBIRHQB3iL1XIhNFkdEcfI6MyfRm5mZNVF5kyxExIKIWFRHsx6kZsB8NxlseS9wYrL9UxFRnrR7CWibu2jNzMy2HXmTLGRpP2Bp2npZUlbVj4En0tYLJb0mabakI3IZoJmZWVOzWbdObilJM0nNYFnVZRHxSDa7qKEsqnzGZUA5/3665DKgfUSsktQdeFhS54j4tIb4TgdOB2jfvn0W4ZiZmTV9DZosRET/LdxFGdAubb0t8MGmFUljgEHAURERyWeuA9Yly3MlLQYOJPVwqarxTQGmAJSUlETVejMzs23R1nYZ4hXgAEmFkr4BjASmQ+ouCeBi4ISIWLtpA0l7JgMjkbQ/cADwboNHbmZmtpXKm2RB0hBJZaSmwn5c0pNJ+b6SZgAkAxjPBp4EFpCaKnt+sosbgVbA01Vukfwu8LqkvwH3A2emPb7azMzM6qCkt96qKCkpiTlzql2pMDMza7IkzY2Ias86ypueBTMzM8tPThbMzMwsIycLZmZmlpGTBTMzM8vIyYKZmZll5GTBzMzMMnKyYGZmZhk5WTAzM7OMnCyYmZlZRk4WzMzMLCMnC2ZmZpaRkwUzMzPLyMmCmZmZZeRkwczMzDJysmBmZmYZOVkwMzOzjPImWZA0TNJ8SRsllWRoN0DSIknvSJqQVn6lpPclzUteA9PqLknaL5J0TK6/i5mZWVPSvLEDSPMGcBLwm9oaSGoG3AR8DygDXpE0PSLeTJpcFxHXVNmmEzAS6AzsC8yUdGBEbMjBdzAzM2ty8qZnISIWRMSiOpr1AN6JiHcj4ivgXuDEOrY5Ebg3ItZFxHvAO8l+zMzMLAt5kyxkaT9gadp6WVK2ydmSXpd0u6TdstzGzMzMMmjQZEHSTElv1PCqq3egYhc1lEXyPhn4FlAMLAOuzWKbqvGdLmmOpDkrVqzIMiQzM7OmrUHHLERE/y3cRRnQLm29LfBBsu8PNxVKuhV4rK5taohvCjAFoKSkpMaEwszMbFuztV2GeAU4QFKhpG+QGrg4HUBSm7R2Q0gNmCSpHylpB0mFwAHAXxswZjMzs61a3iQLkoZIKgN6AY9LejIp31fSDICIKAfOBp4EFgDTImJ+sotfSvq7pNeBfsB5yTbzgWnAm8CfgHG+E8LMzCx7inBve01KSkpizpw5jR2GmZlZg5E0NyKqPevIyUItJK0A/lHPu90DWFnP+7TqfJwbho9zw/Bxbhg+zikdImLPqoVOFhqQpDk1ZWxWv3ycG4aPc8PwcW4YPs6Z5c2YBTMzM8tPThbMzMwsIycLDWtKYwewjfBxbhg+zg3Dx7lh+Dhn4DELZmZmlpF7FszMzCwjJwsNQNIASYskvSNpQmPH0xRJaifpGUkLJM2XdG5jx9SUSWom6TVJj9Xd2r4uSbtKul/SwuTc7tXYMTVFks5L/t94Q9IfJLVo7JjyjZOFHJPUDLgJOBboBHxfUqfGjapJKgcuiIiOQE9gnI9zTp1L6imqllu/Bv4UEQcDXfExr3eS9gP+CyiJiEOAZqSmErA0ThZyrwfwTkS8GxFfAfcC2c6yaVmKiGUR8Wqy/Bmp/1Q9FXkOSGoLHAf8trFjacok7Qx8F7gNICK+iojVjRpU09Uc2FFSc6AltUw2uC1zspB7+wFL09bL8I9YTkkqALoBLzdyKE3V9cB/AxsbOY6mbn9gBXBHcsnnt5K+2dhBNTUR8T5wDfBPYBnwSUQ81bhR5R8nC7mnGsp8C0qOSNoJeAAYHxGfNnY8TY2kQcDyiJjb2LFsA5oDhwKTI6Ib8DngMU/1TNJupHp7C4F9gW9KOqVxo8o/ThZyrwxol7beFndx5YSk7UklClMj4sHGjqeJOhw4QdISUpfUjpR0T+OG1GSVAWURsamH7H5SyYPVr/7AexGxIiLWAw8CvRs5przjZCH3XgEOkFQo6RukBs5Mb+SYmhxJInVtd0FE/Kqx42mqIuKSiGgbEQWkzuU/R4T/CsuBiPgXsFTSQUnRUcCbjRhSU/VPoKeklsn/I0fhgaTVNG/sAJq6iCiXdDbwJKlRtrdHxPxGDqspOhwYDfxd0ryk7NKImNF4IZltsXOAqckfGu8CP2rkeJqciHhZ0v3Aq6TuqnoNP82xGj/B0czMzDLyZQgzMzPLyMmCmZmZZeRkwczMzDJysmBmZmYZOVkwMzOzjJwsmJmZWUZOFszMzCwjJwtm1mgk7ShptqRmkgokvfE19/MNSc8mswaaWT1zsmBmjenHwIMRsWFLdpJM/z4LGFEvUZlZJU4WzCwnJD0j6XvJ8s8k3VBDs1HAIzVsu38yLfNhSY/DwmSK5jckTZXUX9Lzkt6W1CPZ7OFkf2ZWz5wsmFmuXAFcJmkU0A04L70yme9g/4hYUqX8IFKzh/4oIl5Jiv8D+DXQBTgY+AHwHeBC4NKkzRvAYTn5JmbbOF/fM7OciIhnk1n8zgf61nCpYQ9gdZWyPUn1NJxcZcK19yLi7wCS5gOzIiIk/R0oSD5vg6SvJLWKiM/q/xuZbbvcs2BmOSGpCGgDrKvlx/sLoEWVsk+ApaRmEU23Lm15Y9r6Rir/0bMD8OXXjdnMauZkwczqnaQ2wFTgROBzScdUbRMRHwPNJKUnDF8Bg4EfSvrBZn5ma2BFRKz/2oGbWY2cLJhZvZLUEngQuCAiFgBXAVfW0vwpUmMPKkTE58Ag4DxJJ27GR/cDZmx2wGZWJ0VEY8dgZtsoSd2A8yNidD3s60HgkohYtOWRmVk69yyYWaOJiNeAZyQ125L9JHdWPOxEwSw33LNgZmZmGblnwczMzDJysmBmZmYZOVkwMzOzjJwsmJmZWUZOFszMzCwjJwtmZmaWkZMFMzMzy+j/ADAq5QLGdL6rAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x144 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dhdBC3c = np.empty(0)\n",
    "for oc,xp in enumerate(X[::1000]):\n",
    "    M1 = M0\n",
    "    flopy.modflow.ModflowWel(M1, stress_period_data={0: [0, 0, oc, 1.]})\n",
    "    M1.write_input()\n",
    "    success, buff = M1.run_model(silent=True)\n",
    "    A = flopy.utils.binaryfile.HeadFile('model.hds').get_data()[0,0,:]\n",
    "    dhdBC3c = np.append(dhdBC3c, A[-1]*(BC3c-hBC3))\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "plt.subplots(figsize=[8,2])\n",
    "plt.plot(X[::1000]/1000., dhdBC3c, 'k--', mfc='none', label='adjoint sensitivity, numerical')\n",
    "plt.xlabel('$x$ (km)')\n",
    "plt.ylabel(r'$\\partial h / \\partial C_{\\Gamma_3}$')\n",
    "plt.legend();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<img src onerror=\"var cell = this.closest('.code_cell');var editor = cell.querySelector('.input_area');editor.style.background='rgba(0, 200, 0, 0.2)';this.parentNode.removeChild(this)\">"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.000000\n"
     ]
    }
   ],
   "source": [
    "set_background('rgba(0, 200, 0, 0.2)')\n",
    "print('%.6f'% dhdBC3c[int(np.where(X[::1000]==float(ocol))[0])])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}