From fdae4a0e63f5df55de4473302e82509536f1eabe Mon Sep 17 00:00:00 2001
From: Martin Vonk <vonk.mart@gmail.com>
Date: Mon, 9 Jan 2023 18:06:45 +0100
Subject: [PATCH] format all functions and methods according to pep8

---
 pastas/io/men.py           |   8 +-
 pastas/model.py            | 184 +++++++++++++++++++++---------
 pastas/modelcompare.py     |  72 ++++++++----
 pastas/modelplots.py       | 142 ++++++++++++++++-------
 pastas/modelstats.py       | 109 ++++++++++++------
 pastas/noisemodels.py      |   3 +-
 pastas/plots.py            |  55 ++++++---
 pastas/recharge.py         | 116 ++++++++++++++-----
 pastas/reservoir.py        |   7 +-
 pastas/rfunc.py            | 104 ++++++++++++++---
 pastas/solver.py           | 104 ++++++++++++++---
 pastas/stats/core.py       |  40 +++++--
 pastas/stats/dutch.py      |  77 +++++++++----
 pastas/stats/metrics.py    |  63 ++++++++--
 pastas/stats/signatures.py |  49 +++++---
 pastas/stats/tests.py      |  19 +++-
 pastas/stressmodels.py     | 228 ++++++++++++++++++++++++++-----------
 pastas/timer.py            |   5 +-
 pastas/timeseries.py       |   9 +-
 pastas/transform.py        |   8 +-
 pastas/utils.py            |  27 +++--
 21 files changed, 1058 insertions(+), 371 deletions(-)

diff --git a/pastas/io/men.py b/pastas/io/men.py
index 7a6d00628..fc0e6add9 100644
--- a/pastas/io/men.py
+++ b/pastas/io/men.py
@@ -19,8 +19,12 @@ def load(fname: str) -> NotImplementedError:
                               "reads-section for a Menyanthes-read")
 
 
-def dump(fname: str, data: dict, version: int = 3, verbose: bool = True) -> None:
-    # version can also be a specific version, like '2.x.g.t (beta)', or an integer (see below)
+def dump(fname: str,
+         data: dict,
+         version: int = 3,
+         verbose: bool = True) -> None:
+    # version can also be a specific version,
+    # like '2.x.g.t (beta)', or an integer (see below)
     if version == 3:
         version = '3.x.b.c (gamma)'
     elif version == 2:
diff --git a/pastas/model.py b/pastas/model.py
index f8a62cc1d..64de8417b 100644
--- a/pastas/model.py
+++ b/pastas/model.py
@@ -145,7 +145,9 @@ def __repr__(self):
                                const=not self.constant is None,
                                noise=not self.noisemodel is None)
 
-    def add_stressmodel(self, stressmodel: StressModel, replace: bool = False) -> None:
+    def add_stressmodel(self,
+                        stressmodel: StressModel,
+                        replace: bool = False) -> None:
         """Add a stressmodel to the main model.
 
         Parameters
@@ -304,8 +306,12 @@ def del_noisemodel(self) -> None:
             self.noisemodel = None
             self.parameters = self.get_init_parameters(initial=False)
 
-    def simulate(self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None,
-                 tmax: Optional[TimestampType] = None, freq: Optional[str] = None, warmup: Optional[float] = None,
+    def simulate(self,
+                 p: Optional[ArrayLike] = None,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 warmup: Optional[float] = None,
                  return_warmup: bool = False) -> Series:
         """Method to simulate the time series model.
 
@@ -396,9 +402,12 @@ def simulate(self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType]
         sim.name = 'Simulation'
         return sim
 
-    def residuals(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, warmup: Optional[float] = None) -> Series:
+    def residuals(self,
+                  p: Optional[ArrayLike] = None,
+                  tmin: Optional[TimestampType] = None,
+                  tmax: Optional[TimestampType] = None,
+                  freq: Optional[str] = None,
+                  warmup: Optional[float] = None) -> Series:
         """Method to calculate the residual series.
 
         Parameters
@@ -467,9 +476,12 @@ def residuals(
         res.name = "Residuals"
         return res
 
-    def noise(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, warmup: Optional[float] = None) -> Series:
+    def noise(self,
+              p: Optional[ArrayLike] = None,
+              tmin: Optional[TimestampType] = None,
+              tmax: Optional[TimestampType] = None,
+              freq: Optional[str] = None,
+              warmup: Optional[float] = None) -> Series:
         """Method to simulate the noise when a noisemodel is present.
 
         Parameters
@@ -526,9 +538,12 @@ def noise(
         noise = self.noisemodel.simulate(res, p)
         return noise
 
-    def noise_weights(
-            self, p: Optional[list] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, warmup: Optional[float] = None) -> ArrayLike:
+    def noise_weights(self,
+                      p: Optional[list] = None,
+                      tmin: Optional[TimestampType] = None,
+                      tmax: Optional[TimestampType] = None,
+                      freq: Optional[str] = None,
+                      warmup: Optional[float] = None) -> ArrayLike:
         """Internal method to calculate the noise weights."""
         # Get parameters if none are provided
         if p is None:
@@ -542,9 +557,11 @@ def noise_weights(
 
         return weights
 
-    def observations(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, freq: Optional[str] = None,
-            update_observations: bool = False) -> Series:
+    def observations(self,
+                     tmin: Optional[TimestampType] = None,
+                     tmax: Optional[TimestampType] = None,
+                     freq: Optional[str] = None,
+                     update_observations: bool = False) -> Series:
         """Method that returns the observations series used for calibration.
 
         Parameters
@@ -605,10 +622,15 @@ def observations(
             oseries_calib = self.oseries_calib
         return oseries_calib
 
-    def initialize(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, freq: Optional[str] = None,
-            warmup: Optional[float] = None, noise: Optional[bool] = None, weights: Optional[Series] = None, initial:
-            bool = True, fit_constant: bool = True) -> None:
+    def initialize(self,
+                   tmin: Optional[TimestampType] = None,
+                   tmax: Optional[TimestampType] = None,
+                   freq: Optional[str] = None,
+                   warmup: Optional[float] = None,
+                   noise: Optional[bool] = None,
+                   weights: Optional[Series] = None,
+                   initial: bool = True,
+                   fit_constant: bool = True) -> None:
         """Method to initialize the model.
 
         This method is called by the solve-method, but can also be
@@ -662,10 +684,18 @@ def initialize(
             self.parameters.loc["constant_d", "initial"] = 0.0
             self.normalize_residuals = True
 
-    def solve(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, freq: Optional[str] = None,
-            warmup: Optional[float] = None, noise: bool = True, solver: Optional[Solver] = None, report: bool = True,
-            initial: bool = True, weights: Optional[Series] = None, fit_constant: bool = True, **kwargs) -> None:
+    def solve(self,
+              tmin: Optional[TimestampType] = None,
+              tmax: Optional[TimestampType] = None,
+              freq: Optional[str] = None,
+              warmup: Optional[float] = None,
+              noise: bool = True,
+              solver: Optional[Solver] = None,
+              report: bool = True,
+              initial: bool = True,
+              weights: Optional[Series] = None,
+              fit_constant: bool = True,
+              **kwargs) -> None:
         """Method to solve the time series model.
 
         Parameters
@@ -766,9 +796,14 @@ def solve(
                 output = None
             print(self.fit_report(output=output))
 
-    def set_parameter(
-            self, name: str, initial: Optional[float] = None, vary: Optional[bool] = None, pmin: Optional[float] = None,
-            pmax: Optional[float] = None, optimal: Optional[float] = None, move_bounds: bool = False) -> None:
+    def set_parameter(self,
+                      name: str,
+                      initial: Optional[float] = None,
+                      vary: Optional[bool] = None,
+                      pmin: Optional[float] = None,
+                      pmax: Optional[float] = None,
+                      optimal: Optional[float] = None,
+                      move_bounds: bool = False) -> None:
         """Method to change the parameter properties.
 
         Parameters
@@ -881,7 +916,12 @@ def _get_time_offset(self, freq: str) -> Timedelta:
         else:
             return Timedelta(0)
 
-    def _get_sim_index(self, tmin: Timestamp, tmax: Timestamp, freq: str, warmup: Timedelta, update_sim_index: bool = False) -> DatetimeIndex:
+    def _get_sim_index(self,
+                       tmin: Timestamp,
+                       tmax: Timestamp,
+                       freq: str,
+                       warmup: Timedelta,
+                       update_sim_index: bool = False) -> DatetimeIndex:
         """Internal method to get the simulation index, including the warmup.
 
         Parameters
@@ -923,8 +963,10 @@ def _get_sim_index(self, tmin: Timestamp, tmax: Timestamp, freq: str, warmup: Ti
             sim_index = self.sim_index
         return sim_index
 
-    def get_tmin(
-            self, tmin: Optional[TimestampType] = None, use_oseries: bool = True, use_stresses: bool = False) -> Timestamp:
+    def get_tmin(self,
+                 tmin: Optional[TimestampType] = None,
+                 use_oseries: bool = True,
+                 use_stresses: bool = False) -> Timestamp:
         """Method that checks and returns valid values for tmin.
 
         Parameters
@@ -982,8 +1024,10 @@ def get_tmin(
 
         return tmin
 
-    def get_tmax(
-            self, tmax: Optional[TimestampType] = None, use_oseries: bool = True, use_stresses: bool = False) -> Timestamp:
+    def get_tmax(self,
+                 tmax: Optional[TimestampType] = None,
+                 use_oseries: bool = True,
+                 use_stresses: bool = False) -> Timestamp:
         """Method that checks and returns valid values for tmax.
 
         Parameters
@@ -1044,7 +1088,9 @@ def get_tmax(
 
         return tmax
 
-    def get_init_parameters(self, noise: Optional[bool] = None, initial: bool = True) -> DataFrame:
+    def get_init_parameters(self,
+                            noise: Optional[bool] = None,
+                            initial: bool = True) -> DataFrame:
         """Method to get all initial parameters from the individual objects.
 
         Parameters
@@ -1141,10 +1187,15 @@ def get_stressmodel_settings(self, name: str) -> Union[dict, None]:
         return settings
 
     @get_stressmodel
-    def get_contribution(
-            self, name: str, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, warmup: Optional[float] = None, istress: Optional[int] = None,
-            return_warmup: bool = False, p: Optional[ArrayLike] = None) -> Series:
+    def get_contribution(self,
+                         name: str,
+                         tmin: Optional[TimestampType] = None,
+                         tmax: Optional[TimestampType] = None,
+                         freq: Optional[str] = None,
+                         warmup: Optional[float] = None,
+                         istress: Optional[int] = None,
+                         return_warmup: bool = False,
+                         p: Optional[ArrayLike] = None) -> Series:
         """Method to get the contribution of a stressmodel.
 
         Parameters
@@ -1247,7 +1298,9 @@ def get_contributions(self, split: bool = True, **kwargs) -> List[Series]:
         return contribs
 
     def get_transform_contribution(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None) -> Series:
+            self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None) -> Series:
         """Method to get the contribution of a transform.
 
         Parameters
@@ -1271,9 +1324,11 @@ def get_transform_contribution(
         sim_org = ml.simulate(tmin=tmin, tmax=tmax)
         return sim - sim_org
 
-    def get_output_series(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, add_contributions: bool = True,
-            split: bool = True) -> DataFrame:
+    def get_output_series(self,
+                          tmin: Optional[TimestampType] = None,
+                          tmax: Optional[TimestampType] = None,
+                          add_contributions: bool = True,
+                          split: bool = True) -> DataFrame:
         """Method to get all the modeled output time series from the Model.
 
         Parameters
@@ -1324,9 +1379,13 @@ def get_output_series(
         df = concat(df, axis=1)
         return df
 
-    def _get_response(
-            self, block_or_step: str, name: str, p: Optional[ArrayLike] = None, dt: Optional[float] = None, add_0:
-            bool = False, **kwargs) -> Series:
+    def _get_response(self,
+                      block_or_step: str,
+                      name: str,
+                      p: Optional[ArrayLike] = None,
+                      dt: Optional[float] = None,
+                      add_0: bool = False,
+                      **kwargs) -> Series:
         """Internal method to compute the block and step response.
 
         Parameters
@@ -1381,8 +1440,12 @@ def _get_response(
         return response
 
     @get_stressmodel
-    def get_block_response(self, name: str, p: Optional[ArrayLike] = None, add_0: bool = False, dt: Optional[float] = None, **
-                           kwargs) -> Series:
+    def get_block_response(self,
+                           name: str,
+                           p: Optional[ArrayLike] = None,
+                           add_0: bool = False,
+                           dt: Optional[float] = None,
+                           **kwargs) -> Series:
         """Method to obtain the block response for a stressmodel.
 
         The optimal parameters are used when available, initial otherwise.
@@ -1411,8 +1474,12 @@ def get_block_response(self, name: str, p: Optional[ArrayLike] = None, add_0: bo
                                   p=p, add_0=add_0, **kwargs)
 
     @get_stressmodel
-    def get_step_response(
-            self, name, p: Optional[ArrayLike] = None, add_0: bool = False, dt: Optional[float] = None, **kwargs) -> Series:
+    def get_step_response(self,
+                          name,
+                          p: Optional[ArrayLike] = None,
+                          add_0: bool = False,
+                          dt: Optional[float] = None,
+                          **kwargs) -> Series:
         """Method to obtain the step response for a stressmodel.
 
         The optimal parameters are used when available, initial otherwise.
@@ -1440,7 +1507,10 @@ def get_step_response(
                                   p=p, add_0=add_0, **kwargs)
 
     @get_stressmodel
-    def get_response_tmax(self, name: str, p: ArrayLike = None, cutoff: float = 0.999) -> float:
+    def get_response_tmax(self,
+                          name: str,
+                          p: ArrayLike = None,
+                          cutoff: float = 0.999) -> float:
         """Method to get the tmax used for the response function.
 
         Parameters
@@ -1477,9 +1547,16 @@ def get_response_tmax(self, name: str, p: ArrayLike = None, cutoff: float = 0.99
             return tmax
 
     @get_stressmodel
-    def get_stress(self, name: str, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-                   freq: Optional[str] = None, warmup: Optional[float] = None, istress: Optional[int] = None,
-                   return_warmup: bool = False, p: Optional[ArrayLike] = None) -> Union[Series, List[Series]]:
+    def get_stress(
+            self,
+            name: str,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None,
+            freq: Optional[str] = None,
+            warmup: Optional[float] = None,
+            istress: Optional[int] = None,
+            return_warmup: bool = False,
+            p: Optional[ArrayLike] = None) -> Union[Series, List[Series]]:
         """Method to obtain the stress(es) from the stressmodel.
 
         Parameters
@@ -1569,7 +1646,10 @@ def _get_file_info(self) -> dict:
 
         return file_info
 
-    def fit_report(self, output: str = "basic", warnings: bool = True, warnbounds: Optional[bool] = None) -> str:
+    def fit_report(self,
+                   output: str = "basic",
+                   warnings: bool = True,
+                   warnbounds: Optional[bool] = None) -> str:
         """Method that reports on the fit after a model is optimized.
 
         Parameters
diff --git a/pastas/modelcompare.py b/pastas/modelcompare.py
index c1b0a2268..ab50ccb84 100644
--- a/pastas/modelcompare.py
+++ b/pastas/modelcompare.py
@@ -73,7 +73,10 @@ def __init__(self, models: Optional[List[Model]] = None) -> None:
         self.adjust_height = False
         self.smdict = None
 
-    def initialize_figure(self, mosaic: Optional[List[List[str]]] = None, figsize: Tuple[int] = (10, 8), cmap: str = "tab10") -> None:
+    def initialize_figure(self,
+                          mosaic: Optional[List[List[str]]] = None,
+                          figsize: Tuple[int] = (10, 8),
+                          cmap: str = "tab10") -> None:
         """initialize a custom figure based on a mosaic.
 
         Parameters
@@ -94,7 +97,12 @@ def initialize_figure(self, mosaic: Optional[List[List[str]]] = None, figsize: T
         self.axes = axes
         self.cmap = plt.get_cmap(cmap)
 
-    def initialize_adjust_height_figure(self, mosaic: Optional[List[List[str]]] = None, figsize: Tuple[int] = (10, 8), cmap: str = "tab10", smdict: Optional[dict] = None) -> None:
+    def initialize_adjust_height_figure(self,
+                                        mosaic: Optional[List[
+                                            List[str]]] = None,
+                                        figsize: Tuple[int] = (10, 8),
+                                        cmap: str = "tab10",
+                                        smdict: Optional[dict] = None) -> None:
         """initialize subplots based on a mosaic with equal vertical scales.
 
         The height of each subplot is calculated based on the y-data limits in
@@ -233,7 +241,9 @@ def get_unique_stressmodels(self, models: List[Model] = None) -> List[str]:
 
         return sm_unique
 
-    def get_default_mosaic(self, n_stressmodels: Optional[int] = None) -> List[List[str]]:
+    def get_default_mosaic(self,
+                           n_stressmodels: Optional[int] = None
+                          ) -> List[List[str]]:
         """Get default mosaic for matplotlib.subplot_mosaic().
 
         Parameters
@@ -282,7 +292,9 @@ def get_tmin_tmax(self, models: List[Model] = None) -> DataFrame:
 
         return tmintmax
 
-    def get_metrics(self, models: Optional[List[Model]] = None, metric_selection: Optional[List[str]] = None) -> DataFrame:
+    def get_metrics(self,
+                    models: Optional[List[Model]] = None,
+                    metric_selection: Optional[List[str]] = None) -> DataFrame:
         """get metrics of all models in a DataFrame.
 
         Parameters
@@ -311,8 +323,10 @@ def get_metrics(self, models: Optional[List[Model]] = None, metric_selection: Op
         return metrics
 
     def get_parameters(
-        self, models: Optional[List[Model]] = None, param_col: str = "optimal", param_selection: Optional[List[str]] = None
-    ) -> DataFrame:
+            self,
+            models: Optional[List[Model]] = None,
+            param_col: str = "optimal",
+            param_selection: Optional[List[str]] = None) -> DataFrame:
         """get parameter values of all models in a DataFrame.
 
         Parameters
@@ -347,7 +361,9 @@ def get_parameters(
         else:
             return params
 
-    def get_diagnostics(self, models: Optional[List[Model]] = None, diag_col: str = "P-value") -> DataFrame:
+    def get_diagnostics(self,
+                        models: Optional[List[Model]] = None,
+                        diag_col: str = "P-value") -> DataFrame:
         """Get p-values of statistical tests in a DataFrame.
 
         Parameters
@@ -441,7 +457,7 @@ def plot_residuals(self, axn: str = "res") -> None:
             self.axes[axn].plot(
                 residuals.index,
                 residuals.values,
-                label=f"Residuals",
+                label="Residuals",
                 color=self.cmap(i),
             )
 
@@ -462,12 +478,15 @@ def plot_noise(self, axn: str = "res") -> None:
                 self.axes[axn].plot(
                     noise.index,
                     noise.values,
-                    label=f"Noise",
+                    label="Noise",
                     linestyle="--",
                     color=f"C{i}",
                 )
 
-    def plot_response(self, smdict: Optional[dict] = None, axn: str = "rf{i}", response: str = "step") -> None:
+    def plot_response(self,
+                      smdict: Optional[dict] = None,
+                      axn: str = "rf{i}",
+                      response: str = "step") -> None:
         """plot step or block responses.
 
         Parameters
@@ -505,7 +524,7 @@ def plot_response(self, smdict: Optional[dict] = None, axn: str = "rf{i}", respo
             for j, namlist in self.smdict.items():
                 for smn in namlist:
                     # skip if contribution not in model
-                    if not smn in ml.stressmodels:
+                    if smn not in ml.stressmodels:
                         continue
                     if response == "step":
                         step = ml.get_step_response(smn, add_0=True)
@@ -524,7 +543,10 @@ def plot_response(self, smdict: Optional[dict] = None, axn: str = "rf{i}", respo
                             color=self.cmap(i),
                         )
 
-    def plot_contribution(self, smdict: Optional[dict] = None, axn: str = "con{i}", normalized: bool = False) -> None:
+    def plot_contribution(self,
+                          smdict: Optional[dict] = None,
+                          axn: str = "con{i}",
+                          normalized: bool = False) -> None:
         """plot stressmodel contributions.
 
         Parameters
@@ -566,7 +588,7 @@ def plot_contribution(self, smdict: Optional[dict] = None, axn: str = "con{i}",
         for i, ml in enumerate(self.models):
             for j, namlist in self.smdict.items():
                 for smn in namlist:
-                    if not smn in ml.stressmodels:
+                    if smn not in ml.stressmodels:
                         continue
                     for con in ml.get_contributions(split=False):
                         if smn in con.name:
@@ -585,7 +607,9 @@ def plot_contribution(self, smdict: Optional[dict] = None, axn: str = "con{i}",
                                 color=self.cmap(i),
                             )
 
-    def plot_stress(self, axn: str = "stress", names: Optional[List[str]] = None) -> None:
+    def plot_stress(self,
+                    axn: str = "stress",
+                    names: Optional[List[str]] = None) -> None:
         """plot stresses time series.
 
         Parameters
@@ -643,7 +667,9 @@ def plot_acf(self, axn: str = "acf") -> None:
                 label=label,
             )
 
-    def plot_table(self, axn: str = "table", df: Optional[DataFrame] = None) -> None:
+    def plot_table(self,
+                   axn: str = "table",
+                   df: Optional[DataFrame] = None) -> None:
         """plot dataframe as table
 
         Parameters
@@ -670,9 +696,10 @@ def plot_table(self, axn: str = "table", df: Optional[DataFrame] = None) -> None
         self.axes[axn].set_xticks([])
         self.axes[axn].set_yticks([])
 
-    def plot_table_params(
-        self, axn: str = "tab", param_col: str = "optimal", param_selection: Optional[List[str]] = None
-    ) -> None:
+    def plot_table_params(self,
+                          axn: str = "tab",
+                          param_col: str = "optimal",
+                          param_selection: Optional[List[str]] = None) -> None:
         """plot model parameters table.
 
         Parameters
@@ -699,7 +726,10 @@ def plot_table_params(
         cols = params.columns.to_list()[-1:] + params.columns.to_list()[:-1]
         self.plot_table(axn=axn, df=params[cols])
 
-    def plot_table_metrics(self, axn: str = "met", metric_selection: Optional[List[str]] = None) -> None:
+    def plot_table_metrics(
+            self,
+            axn: str = "met",
+            metric_selection: Optional[List[str]] = None) -> None:
         """plot metrics table.
 
         Parameters
@@ -733,7 +763,9 @@ def plot_table_metrics(self, axn: str = "met", metric_selection: Optional[List[s
         cols = metrics.columns.to_list()[-1:] + metrics.columns.to_list()[:-1]
         self.plot_table(axn=axn, df=metrics[cols].round(2))
 
-    def plot_table_diagnostics(self, axn: str = "diag", diag_col: str = "P-value") -> None:
+    def plot_table_diagnostics(self,
+                               axn: str = "diag",
+                               diag_col: str = "P-value") -> None:
         """plot diagnostics table.
 
         Parameters
diff --git a/pastas/modelplots.py b/pastas/modelplots.py
index 89326a544..d1eba0a71 100644
--- a/pastas/modelplots.py
+++ b/pastas/modelplots.py
@@ -43,8 +43,15 @@ def __repr__(self) -> str:
         return msg
 
     @model_tmin_tmax
-    def plot(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, oseries: bool = True, simulation: bool = True,
-             ax: Optional[Axes] = None, figsize: Optional[tuple] = None, legend: bool = True, **kwargs) -> Axes:
+    def plot(self,
+             tmin: Optional[TimestampType] = None,
+             tmax: Optional[TimestampType] = None,
+             oseries: bool = True,
+             simulation: bool = True,
+             ax: Optional[Axes] = None,
+             figsize: Optional[tuple] = None,
+             legend: bool = True,
+             **kwargs) -> Axes:
         """Make a plot of the observed and simulated series.
 
         Parameters
@@ -101,10 +108,16 @@ def plot(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampTyp
         return ax
 
     @model_tmin_tmax
-    def results(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, figsize: tuple = (10, 8),
-            split: bool = False, adjust_height: bool = True, return_warmup: bool = False, block_or_step: str = 'step',
-            fig: Optional[Figure] = None, **kwargs) -> Axes:
+    def results(self,
+                tmin: Optional[TimestampType] = None,
+                tmax: Optional[TimestampType] = None,
+                figsize: tuple = (10, 8),
+                split: bool = False,
+                adjust_height: bool = True,
+                return_warmup: bool = False,
+                block_or_step: str = 'step',
+                fig: Optional[Figure] = None,
+                **kwargs) -> Axes:
         """Plot different results in one window to get a quick overview.
 
         Parameters
@@ -140,13 +153,15 @@ def results(
         o = self.ml.observations(tmin=tmin, tmax=tmax)
         o_nu = self.ml.oseries.series.drop(o.index)
         if return_warmup:
-            o_nu = o_nu[tmin - self.ml.settings['warmup']: tmax]
+            o_nu = o_nu[tmin - self.ml.settings['warmup']:tmax]
         else:
-            o_nu = o_nu[tmin: tmax]
-        sim = self.ml.simulate(tmin=tmin, tmax=tmax,
+            o_nu = o_nu[tmin:tmax]
+        sim = self.ml.simulate(tmin=tmin,
+                               tmax=tmax,
                                return_warmup=return_warmup)
         res = self.ml.residuals(tmin=tmin, tmax=tmax)
-        contribs = self.ml.get_contributions(split=split, tmin=tmin,
+        contribs = self.ml.get_contributions(split=split,
+                                             tmin=tmin,
                                              tmax=tmax,
                                              return_warmup=return_warmup)
 
@@ -283,11 +298,17 @@ def results(
         return fig.axes
 
     @model_tmin_tmax
-    def decomposition(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, ytick_base: bool = True,
-            split: bool = True, figsize: tuple = (10, 8),
-            axes: Optional[Axes] = None, name: Optional[str] = None, return_warmup: bool = False,
-            min_ylim_diff: Optional[float] = None, **kwargs) -> Axes:
+    def decomposition(self,
+                      tmin: Optional[TimestampType] = None,
+                      tmax: Optional[TimestampType] = None,
+                      ytick_base: bool = True,
+                      split: bool = True,
+                      figsize: tuple = (10, 8),
+                      axes: Optional[Axes] = None,
+                      name: Optional[str] = None,
+                      return_warmup: bool = False,
+                      min_ylim_diff: Optional[float] = None,
+                      **kwargs) -> Axes:
         """Plot the decomposition of a time-series in the different stresses.
 
         Parameters
@@ -419,9 +440,15 @@ def decomposition(
         return axes
 
     @model_tmin_tmax
-    def diagnostics(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-                    figsize: tuple = (10, 5), bins: int = 50, acf_options: Optional[dict] = None,
-                    fig: Optional[Figure] = None, alpha: float = 0.05, **kwargs) -> Axes:
+    def diagnostics(self,
+                    tmin: Optional[TimestampType] = None,
+                    tmax: Optional[TimestampType] = None,
+                    figsize: tuple = (10, 5),
+                    bins: int = 50,
+                    acf_options: Optional[dict] = None,
+                    fig: Optional[Figure] = None,
+                    alpha: float = 0.05,
+                    **kwargs) -> Axes:
         """Plot a window that helps in diagnosing basic model assumptions.
 
         Parameters
@@ -482,10 +509,12 @@ def diagnostics(self, tmin: Optional[TimestampType] = None, tmax: Optional[Times
                            **kwargs)
 
     @model_tmin_tmax
-    def cum_frequency(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, ax: Optional[Axes] = None,
-            figsize: tuple = (5, 2),
-            **kwargs) -> Axes:
+    def cum_frequency(self,
+                      tmin: Optional[TimestampType] = None,
+                      tmax: Optional[TimestampType] = None,
+                      ax: Optional[Axes] = None,
+                      figsize: tuple = (5, 2),
+                      **kwargs) -> Axes:
         """Plot the cumulative frequency for the observations and simulation.
 
         Parameters
@@ -513,9 +542,11 @@ def cum_frequency(
         obs = self.ml.observations(tmin=tmin, tmax=tmax)
         return cum_frequency(obs, sim, ax=ax, figsize=figsize, **kwargs)
 
-    def block_response(
-            self, stressmodels: Optional[List[str]] = None, ax: Optional[Axes] = None, figsize: Optional[tuple] = None, **
-            kwargs) -> Axes:
+    def block_response(self,
+                       stressmodels: Optional[List[str]] = None,
+                       ax: Optional[Axes] = None,
+                       figsize: Optional[tuple] = None,
+                       **kwargs) -> Axes:
         """Plot the block response for a specific stressmodels.
 
         Parameters
@@ -553,9 +584,11 @@ def block_response(
         plt.legend(legend)
         return ax
 
-    def step_response(
-            self, stressmodels: Optional[List[str]] = None, ax: Optional[Axes] = None, figsize: Optional[tuple] = None, **
-            kwargs) -> Axes:
+    def step_response(self,
+                      stressmodels: Optional[List[str]] = None,
+                      ax: Optional[Axes] = None,
+                      figsize: Optional[tuple] = None,
+                      **kwargs) -> Axes:
         """Plot the step response for a specific stressmodels.
 
         Parameters
@@ -590,10 +623,14 @@ def step_response(
         return ax
 
     @model_tmin_tmax
-    def stresses(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, cols: int = 1, split: bool = True,
-            sharex: bool = True, figsize: tuple = (10, 8),
-            **kwargs) -> Axes:
+    def stresses(self,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 cols: int = 1,
+                 split: bool = True,
+                 sharex: bool = True,
+                 figsize: tuple = (10, 8),
+                 **kwargs) -> Axes:
         """This method creates a graph with all the stresses used in the model.
 
         Parameters
@@ -637,10 +674,17 @@ def stresses(
         return axes
 
     @model_tmin_tmax
-    def contributions_pie(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, ax: Optional[Axes] = None,
-            figsize: Optional[Figure] = None, split: bool = True, partition: str = 'std', wedgeprops: Optional[dict] = None,
-            startangle: float = 90.0, autopct: str = "%1.1f%%", **kwargs) -> Axes:
+    def contributions_pie(self,
+                          tmin: Optional[TimestampType] = None,
+                          tmax: Optional[TimestampType] = None,
+                          ax: Optional[Axes] = None,
+                          figsize: Optional[Figure] = None,
+                          split: bool = True,
+                          partition: str = 'std',
+                          wedgeprops: Optional[dict] = None,
+                          startangle: float = 90.0,
+                          autopct: str = "%1.1f%%",
+                          **kwargs) -> Axes:
         """Make a pie chart of the contributions. This plot is based on the TNO
         Groundwatertoolbox.
 
@@ -705,9 +749,13 @@ def contributions_pie(
         return ax
 
     @model_tmin_tmax
-    def stacked_results(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, figsize: tuple = (10, 8),
-            stacklegend: bool = False, stacklegend_kws: Optional[dict] = None, **kwargs) -> Axes:
+    def stacked_results(self,
+                        tmin: Optional[TimestampType] = None,
+                        tmax: Optional[TimestampType] = None,
+                        figsize: tuple = (10, 8),
+                        stacklegend: bool = False,
+                        stacklegend_kws: Optional[dict] = None,
+                        **kwargs) -> Axes:
         """Create a results plot, similar to `ml.plots.results()`, in which the
         individual contributions of stresses (in stressmodels with multiple
         stresses) are stacked.
@@ -789,7 +837,11 @@ def custom_sort(t):
         return axes
 
     @model_tmin_tmax
-    def series(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, split: bool = True, **kwargs) -> Axes:
+    def series(self,
+               tmin: Optional[TimestampType] = None,
+               tmax: Optional[TimestampType] = None,
+               split: bool = True,
+               **kwargs) -> Axes:
         """Method to plot all the time series going into a Pastas Model.
 
         Parameters
@@ -821,9 +873,13 @@ def series(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampT
         return axes
 
     @model_tmin_tmax
-    def summary_pdf(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, fname: Optional[str] = None,
-            dpi: int = 150, results_kwargs: Optional[dict] = None, diagnostics_kwargs: Optional[dict] = None) -> Figure:
+    def summary_pdf(self,
+                    tmin: Optional[TimestampType] = None,
+                    tmax: Optional[TimestampType] = None,
+                    fname: Optional[str] = None,
+                    dpi: int = 150,
+                    results_kwargs: Optional[dict] = None,
+                    diagnostics_kwargs: Optional[dict] = None) -> Figure:
         """Create a PDF file (A4) with the results and diagnostics plot.
 
         Parameters
diff --git a/pastas/modelstats.py b/pastas/modelstats.py
index 8dcba690a..d15e34260 100644
--- a/pastas/modelstats.py
+++ b/pastas/modelstats.py
@@ -35,7 +35,17 @@
 
 class Statistics:
     # Save all statistics that can be calculated.
-    ops = ["rmse", "rmsn", "sse", "mae", "nse", "evp", "rsq", "bic", "aic", ]
+    ops = [
+        "rmse",
+        "rmsn",
+        "sse",
+        "mae",
+        "nse",
+        "evp",
+        "rsq",
+        "bic",
+        "aic",
+    ]
 
     def __init__(self, ml: Model):
         """This class provides statistics to to pastas Model class.
@@ -62,9 +72,11 @@ def __repr__(self):
         return msg
 
     @model_tmin_tmax
-    def rmse(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def rmse(self,
+             tmin: Optional[TimestampType] = None,
+             tmax: Optional[TimestampType] = None,
+             weighted: bool = False,
+             **kwargs) -> float:
         """Root mean squared error of the residuals.
 
         Parameters
@@ -83,9 +95,11 @@ def rmse(
         return metrics.rmse(res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def rmsn(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def rmsn(self,
+             tmin: Optional[TimestampType] = None,
+             tmax: Optional[TimestampType] = None,
+             weighted: bool = False,
+             **kwargs) -> float:
         """Root mean squared error of the noise.
 
         Parameters
@@ -112,7 +126,9 @@ def rmsn(
             return metrics.rmse(res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def sse(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None) -> float:
+    def sse(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None) -> float:
         """Sum of the squares of the error (SSE)
 
         Parameters
@@ -128,9 +144,11 @@ def sse(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType
         return metrics.sse(res=res)
 
     @model_tmin_tmax
-    def mae(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def mae(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None,
+            weighted: bool = False,
+            **kwargs) -> float:
         """Mean Absolute Error (MAE) of the residuals.
 
         Parameters
@@ -149,9 +167,11 @@ def mae(
         return metrics.mae(res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def nse(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def nse(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None,
+            weighted: bool = False,
+            **kwargs) -> float:
         """Nash-Sutcliffe coefficient for model fit .
 
         Parameters
@@ -171,9 +191,11 @@ def nse(
         return metrics.nse(obs=obs, res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def pearsonr(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def pearsonr(self,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 weighted: bool = False,
+                 **kwargs) -> float:
         """Compute the (weighted) Pearson correlation (r).
 
         Parameters
@@ -193,9 +215,11 @@ def pearsonr(
         return metrics.pearsonr(obs=obs, sim=sim, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def evp(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def evp(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None,
+            weighted: bool = False,
+            **kwargs) -> float:
         """Explained variance percentage.
 
         Parameters
@@ -215,9 +239,11 @@ def evp(
         return metrics.evp(obs=obs, res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def rsq(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def rsq(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None,
+            weighted: bool = False,
+            **kwargs) -> float:
         """R-squared.
 
         Parameters
@@ -237,9 +263,11 @@ def rsq(
         return metrics.rsq(obs=obs, res=res, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def kge_2012(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, weighted: bool = False, **
-            kwargs) -> float:
+    def kge_2012(self,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 weighted: bool = False,
+                 **kwargs) -> float:
         """Kling-Gupta Efficiency.
 
         Parameters
@@ -259,7 +287,9 @@ def kge_2012(
         return metrics.kge_2012(obs=obs, sim=sim, weighted=weighted, **kwargs)
 
     @model_tmin_tmax
-    def bic(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None) -> float:
+    def bic(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None) -> float:
         """Bayesian Information Criterium (BIC).
 
         Parameters
@@ -281,7 +311,9 @@ def bic(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType
         return metrics.bic(res=res, nparam=nparam)
 
     @model_tmin_tmax
-    def aic(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None) -> float:
+    def aic(self,
+            tmin: Optional[TimestampType] = None,
+            tmax: Optional[TimestampType] = None) -> float:
         """Akaike Information Criterium (AIC).
 
         Parameters
@@ -302,7 +334,9 @@ def aic(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType
         return metrics.aic(res=res, nparam=nparam)
 
     @model_tmin_tmax
-    def summary(self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
+    def summary(self,
+                tmin: Optional[TimestampType] = None,
+                tmax: Optional[TimestampType] = None,
                 stats: Optional[List[str]] = None) -> DataFrame:
         """Returns a Pandas DataFrame with goodness-of-fit metrics.
 
@@ -342,10 +376,12 @@ def summary(self, tmin: Optional[TimestampType] = None, tmax: Optional[Timestamp
         return stats
 
     @model_tmin_tmax
-    def diagnostics(
-            self, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, alpha: float = 0.05,
-            stats: tuple = (),
-            float_fmt: str = "{0:.2f}") -> DataFrame:
+    def diagnostics(self,
+                    tmin: Optional[TimestampType] = None,
+                    tmax: Optional[TimestampType] = None,
+                    alpha: float = 0.05,
+                    stats: tuple = (),
+                    float_fmt: str = "{0:.2f}") -> DataFrame:
         if self.ml.noisemodel and self.ml.settings["noise"]:
             series = self.ml.noise(tmin=tmin, tmax=tmax)
             nparam = self.ml.noisemodel.nparam
@@ -353,5 +389,8 @@ def diagnostics(
             series = self.ml.residuals(tmin=tmin, tmax=tmax)
             nparam = 0
 
-        return diagnostics(series=series, alpha=alpha, nparam=nparam,
-                           stats=stats, float_fmt=float_fmt)
+        return diagnostics(series=series,
+                           alpha=alpha,
+                           nparam=nparam,
+                           stats=stats,
+                           float_fmt=float_fmt)
diff --git a/pastas/noisemodels.py b/pastas/noisemodels.py
index 57c2e8713..eeb1374da 100644
--- a/pastas/noisemodels.py
+++ b/pastas/noisemodels.py
@@ -246,7 +246,8 @@ def simulate(self, res: Series, p: ArrayLike) -> Series:
 
     @staticmethod
     @njit
-    def calculate_noise(res: ArrayLike, odelt: ArrayLike, alpha: float, beta: float) -> ArrayLike:
+    def calculate_noise(res: ArrayLike, odelt: ArrayLike, alpha: float,
+                        beta: float) -> ArrayLike:
         # Create an array to store the noise
         a = np.zeros_like(res)
         a[0] = res[0]
diff --git a/pastas/plots.py b/pastas/plots.py
index a1ae26887..fd464dff0 100644
--- a/pastas/plots.py
+++ b/pastas/plots.py
@@ -58,9 +58,15 @@ def compare(models: List[Model], adjust_height: bool = True, **kwargs) -> Axes:
 
 
 def series(
-        head: Optional[Series] = None, stresses: Optional[List[Series]] = None, hist: bool = True, kde: bool = False,
-        titles: bool = True, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-        labels: Optional[List[str]] = None, figsize: tuple = (10, 5)) -> Axes:
+    head: Optional[Series] = None,
+    stresses: Optional[List[Series]] = None,
+    hist: bool = True,
+    kde: bool = False,
+    titles: bool = True,
+    tmin: Optional[TimestampType] = None,
+    tmax: Optional[TimestampType] = None,
+    labels: Optional[List[str]] = None,
+    figsize: tuple = (10, 5)) -> Axes:
     """Plot all the input time Series in a single plot.
 
     Parameters
@@ -200,9 +206,14 @@ def series(
     return axes
 
 
-def acf(
-        series: Series, alpha: float = 0.05, lags: int = 365, acf_options: Optional[dict] = None, smooth_conf: bool = True,
-        color: str = "k", ax: Optional[Axes] = None, figsize: tuple = (5, 2)) -> Axes:
+def acf(series: Series,
+        alpha: float = 0.05,
+        lags: int = 365,
+        acf_options: Optional[dict] = None,
+        smooth_conf: bool = True,
+        color: str = "k",
+        ax: Optional[Axes] = None,
+        figsize: tuple = (5, 2)) -> Axes:
     """Plot of the autocorrelation function of a time series.
 
     Parameters
@@ -268,10 +279,15 @@ def acf(
     return ax
 
 
-def diagnostics(
-        series: Series, sim: Optional[Series] = None, alpha: float = 0.05, bins: int = 50, acf_options: Optional[dict] = None,
-        figsize: tuple = (10, 5),
-        fig: Optional[Figure] = None, heteroscedasicity: bool = True, **kwargs) -> Axes:
+def diagnostics(series: Series,
+                sim: Optional[Series] = None,
+                alpha: float = 0.05,
+                bins: int = 50,
+                acf_options: Optional[dict] = None,
+                figsize: tuple = (10, 5),
+                fig: Optional[Figure] = None,
+                heteroscedasicity: bool = True,
+                **kwargs) -> Axes:
     """Plot that helps in diagnosing basic model assumptions.
 
     Parameters
@@ -391,8 +407,10 @@ def diagnostics(
     return fig.axes
 
 
-def cum_frequency(
-        obs: Series, sim: Optional[Series] = None, ax: Optional[Axes] = None, figsize: tuple = (5, 2)) -> Axes:
+def cum_frequency(obs: Series,
+                  sim: Optional[Series] = None,
+                  ax: Optional[Axes] = None,
+                  figsize: tuple = (5, 2)) -> Axes:
     """Plot of the cumulative frequency of a time Series.
 
     Parameters
@@ -495,9 +513,12 @@ class TrackSolve:
     Access the resulting figure through `track.fig`.
     """
 
-    def __init__(
-            self, ml: Model, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            update_iter: Optional[int] = None) -> None:
+    def __init__(self,
+                 ml: Model,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 update_iter: Optional[int] = None) -> None:
+
         logger.warning("TrackSolve feature under development. If you find any "
                        "bugs please post an issue on GitHub: "
                        "https://github.com/pastas/pastas/issues")
@@ -637,7 +658,9 @@ def _simulate(self) -> Series:
                                freq=self.ml.settings["freq"])
         return sim
 
-    def initialize_figure(self, figsize: Tuple[int] = (10, 8), dpi: int = 100) -> Figure:
+    def initialize_figure(self,
+                          figsize: Tuple[int] = (10, 8),
+                          dpi: int = 100) -> Figure:
         """Initialize figure for plotting optimization progress.
 
         Parameters
diff --git a/pastas/recharge.py b/pastas/recharge.py
index 7687e59c8..147b92247 100644
--- a/pastas/recharge.py
+++ b/pastas/recharge.py
@@ -108,7 +108,8 @@ def get_init_parameters(self, name: str = "recharge") -> DataFrame:
         parameters.loc[name + "_f"] = (-1.0, -2.0, 0.0, True, name)
         return parameters
 
-    def simulate(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, **kwargs) -> ArrayLike:
+    def simulate(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike,
+                 **kwargs) -> ArrayLike:
         """Simulate the precipitation excess flux.
 
         Parameters
@@ -128,7 +129,8 @@ def simulate(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, **kwargs) ->
         """
         return add(prec, multiply(evap, p))
 
-    def get_water_balance(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, **kwargs) -> DataFrame:
+    def get_water_balance(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike,
+                          **kwargs) -> DataFrame:
         ea = multiply(evap, p)
         r = add(prec, multiply(evap, p))
         return DataFrame(data=vstack((prec, ea, -r)).T,
@@ -184,7 +186,10 @@ class FlexModel(RechargeBase):
     """
     _name = "FlexModel"
 
-    def __init__(self, interception: bool = True, snow: bool = False, gw_uptake: bool = False):
+    def __init__(self,
+                 interception: bool = True,
+                 snow: bool = False,
+                 gw_uptake: bool = False):
         check_numba()
         RechargeBase.__init__(self)
         self.snow = snow
@@ -216,8 +221,14 @@ def get_init_parameters(self, name: str = "recharge") -> DataFrame:
 
         return parameters
 
-    def simulate(self, prec: ArrayLike, evap: ArrayLike, temp: ArrayLike, p: ArrayLike, dt: float = 1.0,
-                 return_full: bool = False, **kwargs) -> ArrayLike:
+    def simulate(self,
+                 prec: ArrayLike,
+                 evap: ArrayLike,
+                 temp: ArrayLike,
+                 p: ArrayLike,
+                 dt: float = 1.0,
+                 return_full: bool = False,
+                 **kwargs) -> ArrayLike:
         """Simulate the soil water balance model.
 
         Parameters
@@ -295,8 +306,13 @@ def simulate(self, prec: ArrayLike, evap: ArrayLike, temp: ArrayLike, p: ArrayLi
 
     @staticmethod
     @njit
-    def get_root_zone_balance(pe: ArrayLike, ep: ArrayLike, srmax: float = 250.0, lp: float = 0.25, ks: float = 100.0,
-                              gamma: float = 4.0, dt: float = 1.0) -> Tuple[ArrayLike]:
+    def get_root_zone_balance(pe: ArrayLike,
+                              ep: ArrayLike,
+                              srmax: float = 250.0,
+                              lp: float = 0.25,
+                              ks: float = 100.0,
+                              gamma: float = 4.0,
+                              dt: float = 1.0) -> Tuple[ArrayLike]:
         """Method to compute the water balance of the root zone reservoir.
 
         Parameters
@@ -368,8 +384,10 @@ def get_root_zone_balance(pe: ArrayLike, ep: ArrayLike, srmax: float = 250.0, lp
 
     @staticmethod
     @njit
-    def get_interception_balance(
-            pr: ArrayLike, ep: ArrayLike, simax: float = 2.0, dt: float = 1.0) -> Tuple[ArrayLike]:
+    def get_interception_balance(pr: ArrayLike,
+                                 ep: ArrayLike,
+                                 simax: float = 2.0,
+                                 dt: float = 1.0) -> Tuple[ArrayLike]:
         """Method to compute the water balance of the interception reservoir.
 
         Parameters
@@ -423,7 +441,10 @@ def get_interception_balance(
 
     @staticmethod
     @njit
-    def get_snow_balance(prec: ArrayLike, temp: ArrayLike, tt: float = 0.0, k: float = 2.0) -> Tuple[ArrayLike]:
+    def get_snow_balance(prec: ArrayLike,
+                         temp: ArrayLike,
+                         tt: float = 0.0,
+                         k: float = 2.0) -> Tuple[ArrayLike]:
         """Method to compute the water balance of the snow reservoir.
 
         Parameters
@@ -471,8 +492,13 @@ def get_snow_balance(prec: ArrayLike, temp: ArrayLike, tt: float = 0.0, k: float
 
         return ss[:-1], ps, -m
 
-    def get_water_balance(self, prec: ArrayLike, evap: ArrayLike, temp: ArrayLike, p: ArrayLike, dt: float = 1.0, **
-                          kwargs) -> DataFrame:
+    def get_water_balance(self,
+                          prec: ArrayLike,
+                          evap: ArrayLike,
+                          temp: ArrayLike,
+                          p: ArrayLike,
+                          dt: float = 1.0,
+                          **kwargs) -> DataFrame:
         data = self.simulate(prec=prec, evap=evap, temp=temp,
                              p=p, dt=dt, return_full=True, **kwargs)
 
@@ -490,17 +516,20 @@ def get_water_balance(self, prec: ArrayLike, evap: ArrayLike, temp: ArrayLike, p
 
         return DataFrame(data=vstack(data).T, columns=columns)
 
-    def check_snow_balance(self, prec: ArrayLike, temp: ArrayLike, **kwargs) -> float:
+    def check_snow_balance(self, prec: ArrayLike, temp: ArrayLike,
+                           **kwargs) -> float:
         ss, ps, m = self.get_snow_balance(prec, temp)
         error = (ss[0] - ss[-1] + (ps + m).sum())
         return error
 
-    def check_interception_balance(self, prec: ArrayLike, evap: ArrayLike, **kwargs) -> float:
+    def check_interception_balance(self, prec: ArrayLike, evap: ArrayLike,
+                                   **kwargs) -> float:
         si, ei, pi = self.get_interception_balance(prec, evap)
         error = (si[0] - si[-1] + (pi + ei).sum())
         return error
 
-    def check_root_zone_balance(self, prec: ArrayLike, evap: ArrayLike, **kwargs) -> float:
+    def check_root_zone_balance(self, prec: ArrayLike, evap: ArrayLike,
+                                **kwargs) -> float:
         sr, r, ea, q, pe = self.get_root_zone_balance(prec, evap)
         error = (sr[0] - sr[-1] + (r + ea + q + pe).sum())
         return error
@@ -546,9 +575,13 @@ def get_init_parameters(self, name: str = "recharge") -> DataFrame:
         parameters.loc[name + "_ks"] = (100.0, 1, 1e4, True, name)
         return parameters
 
-    def simulate(
-            self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, dt: ArrayLike = 1.0, return_full: bool = False, **
-            kwargs) -> Tuple[ArrayLike]:
+    def simulate(self,
+                 prec: ArrayLike,
+                 evap: ArrayLike,
+                 p: ArrayLike,
+                 dt: ArrayLike = 1.0,
+                 return_full: bool = False,
+                 **kwargs) -> Tuple[ArrayLike]:
         """Simulate the recharge flux.
 
         Parameters
@@ -580,8 +613,16 @@ def simulate(
 
     @staticmethod
     @njit
-    def get_recharge(prec: ArrayLike, evap: ArrayLike, fi: float = 1.0, fc: float = 1.0, sr: float = 0.5, de: float = 250.0,
-                     l: float = -2.0, m: float = 0.5, ks: float = 50.0, dt: float = 1.0) -> ArrayLike:
+    def get_recharge(prec: ArrayLike,
+                     evap: ArrayLike,
+                     fi: float = 1.0,
+                     fc: float = 1.0,
+                     sr: float = 0.5,
+                     de: float = 250.0,
+                     l: float = -2.0,
+                     m: float = 0.5,
+                     ks: float = 50.0,
+                     dt: float = 1.0) -> ArrayLike:
         """
         Internal method used for the recharge calculation. If Numba is
         available, this method is significantly faster.
@@ -613,8 +654,12 @@ def get_recharge(prec: ArrayLike, evap: ArrayLike, fi: float = 1.0, fc: float =
             s[t + 1] = s[t] + dt / de * (pe[t] - ea[t] - r[t])
         return r, s, ea, pe
 
-    def get_water_balance(
-            self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, dt: float = 1.0, **kwargs) -> DataFrame:
+    def get_water_balance(self,
+                          prec: ArrayLike,
+                          evap: ArrayLike,
+                          p: ArrayLike,
+                          dt: float = 1.0,
+                          **kwargs) -> DataFrame:
         r, s, ea, pe = self.simulate(prec, evap, p=p, dt=dt,
                                      return_full=True, **kwargs)
         s = s * p[3]  # Because S is computed dimensionless in this model
@@ -674,8 +719,13 @@ def get_init_parameters(self, name: str = "recharge") -> DataFrame:
         parameters.loc[name + "_gamma"] = (1.0, 0.0, 2.0, True, name)
         return parameters
 
-    def simulate(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, dt: float = 1.0, return_full: bool = False, **
-                 kwargs) -> ArrayLike:
+    def simulate(self,
+                 prec: ArrayLike,
+                 evap: ArrayLike,
+                 p: ArrayLike,
+                 dt: float = 1.0,
+                 return_full: bool = False,
+                 **kwargs) -> ArrayLike:
         """Simulate the recharge flux.
 
         Parameters
@@ -706,8 +756,14 @@ def simulate(self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, dt: float = 1
 
     @staticmethod
     @njit
-    def get_recharge(prec: ArrayLike, evap: ArrayLike, scap: float = 1.0, alpha: float = 1.0,
-                     ksat: float = 1.0, beta: float = 0.5, gamma: float = 1.0, dt: float = 1.0):
+    def get_recharge(prec: ArrayLike,
+                     evap: ArrayLike,
+                     scap: float = 1.0,
+                     alpha: float = 1.0,
+                     ksat: float = 1.0,
+                     beta: float = 0.5,
+                     gamma: float = 1.0,
+                     dt: float = 1.0):
         """
         Internal method used for the recharge calculation. If Numba is
         available, this method is significantly faster.
@@ -735,8 +791,12 @@ def get_recharge(prec: ArrayLike, evap: ArrayLike, scap: float = 1.0, alpha: flo
                             max(0.0, sm[t] + (pe[t] - ea[t] - r[t]) * dt))
         return r, sm[1:], ea, pe
 
-    def get_water_balance(
-            self, prec: ArrayLike, evap: ArrayLike, p: ArrayLike, dt: float = 1.0, **kwargs) -> DataFrame:
+    def get_water_balance(self,
+                          prec: ArrayLike,
+                          evap: ArrayLike,
+                          p: ArrayLike,
+                          dt: float = 1.0,
+                          **kwargs) -> DataFrame:
         r, s, ea, pe = self.simulate(prec, evap, p=p, dt=dt,
                                      return_full=True, **kwargs)
         data = DataFrame(data=vstack((s, pe, ea, r)).T,
diff --git a/pastas/reservoir.py b/pastas/reservoir.py
index 1d6a2932e..2c8ebe153 100644
--- a/pastas/reservoir.py
+++ b/pastas/reservoir.py
@@ -47,7 +47,12 @@ def get_init_parameters(name: str = "recharge") -> DataFrame:
             columns=["initial", "pmin", "pmax", "vary", "name"])
         return parameters
 
-    def simulate(self, prec: Series, evap: Series, p: ArrayLike, dt: float = 1.0, **kwargs) -> ArrayLike:
+    def simulate(self,
+                 prec: Series,
+                 evap: Series,
+                 p: ArrayLike,
+                 dt: float = 1.0,
+                 **kwargs) -> ArrayLike:
         pass
 
 
diff --git a/pastas/rfunc.py b/pastas/rfunc.py
index e9bd018ce..a0f829847 100644
--- a/pastas/rfunc.py
+++ b/pastas/rfunc.py
@@ -37,7 +37,9 @@ def __init__(self, **kwargs) -> None:
         self.cutoff = 0.999
         self.kwargs = kwargs
 
-    def _set_init_parameter_settings(self, up: bool = True, meanstress: float = 1.0,
+    def _set_init_parameter_settings(self,
+                                     up: bool = True,
+                                     meanstress: float = 1.0,
                                      cutoff: float = 0.999) -> None:
         self.up = up
         # Completely arbitrary number to prevent division by zero
@@ -82,7 +84,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
         """
         pass
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         """Method to return the step function.
 
         Parameters
@@ -104,7 +110,11 @@ def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, ma
         """
         pass
 
-    def block(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def block(self,
+              p: ArrayLike,
+              dt: float = 1.0,
+              cutoff: Optional[float] = None,
+              maxtmax: Optional[int] = None) -> ArrayLike:
         """Method to return the block function.
 
         Parameters
@@ -153,7 +163,11 @@ def impulse(self, t: ArrayLike, p: ArrayLike) -> ArrayLike:
         """
         pass
 
-    def get_t(self, p: ArrayLike, dt: float, cutoff: float, maxtmax: Optional[int] = None) -> ArrayLike:
+    def get_t(self,
+              p: ArrayLike,
+              dt: float,
+              cutoff: float,
+              maxtmax: Optional[int] = None) -> ArrayLike:
         """Internal method to determine the times at which to evaluate the
         step-response, from t=0.
 
@@ -242,7 +256,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(self, p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         s = p[0] * gammainc(p[1], t / p[2])
         return s
@@ -306,7 +324,11 @@ def get_tmax(self, p: ArrayLike, cutoff=None) -> float:
     def gain(self, p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[float] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[float] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         s = p[0] * (1.0 - np.exp(-t / p[1]))
         return s
@@ -461,7 +483,8 @@ def numpy_step(A: float, a: float, b: float, r: float, t: ArrayLike) -> ArrayLik
             tau2 + rhosq / (4 * tau2))
         return A * F / 2
 
-    def quad_step(self, A: float, a: float, b: float, r: float, t: ArrayLike) -> ArrayLike:
+    def quad_step(self, A: float, a: float, b: float, r: float,
+                  t: ArrayLike) -> ArrayLike:
         F = np.zeros_like(t)
         brsq = np.exp(b) * r**2
         u = a * brsq / t
@@ -470,7 +493,11 @@ def quad_step(self, A: float, a: float, b: float, r: float, t: ArrayLike) -> Arr
                         u[i], np.inf, args=(brsq,))[0]
         return F * A / 2
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         A, a, b = p[:3]
         r = self._get_distance_from_params(p)
         t = self.get_t(p, dt, cutoff, maxtmax)
@@ -485,7 +512,12 @@ def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, ma
                 return self.numpy_step(A, a, b, r, t)
 
     @staticmethod
-    def variance_gain(A: float, b: float, var_A: float, var_b: float, cov_Ab: float, r: Optional[float] = 1.0) -> Union[float, ArrayLike]:
+    def variance_gain(A: float,
+                      b: float,
+                      var_A: float,
+                      var_b: float,
+                      cov_Ab: float,
+                      r: Optional[float] = 1.0) -> Union[float, ArrayLike]:
         """Calculate variance of the gain from parameters A and b.
 
         Variance of the gain is calculated based on propagation of
@@ -658,7 +690,11 @@ def quad_step(self, A: float, a: float, b: float, t: ArrayLike) -> ArrayLike:
                         u[i], np.inf, args=(b,))[0]
         return F * A / (2 * k0(2 * np.sqrt(b)))
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         A, a, b = p
         t = self.get_t(p, dt, cutoff, maxtmax)
 
@@ -727,7 +763,11 @@ def gain(self, p: ArrayLike) -> float:
             g = -g
         return g
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         A, a, b = p
         s = A * self.polder_function(np.sqrt(b), np.sqrt(t / a))
@@ -788,13 +828,21 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(self, p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         if isinstance(dt, np.ndarray):
             return p[0] * np.ones(len(dt))
         else:
             return p[0] * np.ones(1)
 
-    def block(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def block(self,
+              p: ArrayLike,
+              dt: float = 1.0,
+              cutoff: Optional[float] = None,
+              maxtmax: Optional[int] = None) -> ArrayLike:
         return p[0] * np.ones(1)
 
 
@@ -890,7 +938,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
 
         if self.quad:
             t = self.get_t(p, dt, cutoff, maxtmax)
@@ -1006,7 +1058,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(self, p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         s = p[0] * (1 - ((1 - p[1]) * np.exp(-t / p[2]) +
                          p[1] * np.exp(-t / p[3])))
@@ -1062,7 +1118,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(p: ArrayLike) -> float:
         return 1.
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         s = erfc(1 / (p[0] * np.sqrt(t)))
         return s
@@ -1137,7 +1197,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         t = self.get_t(p, dt, cutoff, maxtmax)
         h = 0
         for n in range(self.n_terms):
@@ -1218,7 +1282,11 @@ def get_tmax(self, p: ArrayLike, cutoff: Optional[float] = None) -> float:
     def gain(self, p: ArrayLike) -> float:
         return p[0]
 
-    def step(self, p: ArrayLike, dt: float = 1.0, cutoff: Optional[float] = None, maxtmax: Optional[int] = None) -> ArrayLike:
+    def step(self,
+             p: ArrayLike,
+             dt: float = 1.0,
+             cutoff: Optional[float] = None,
+             maxtmax: Optional[int] = None) -> ArrayLike:
         f = interp1d(self.t, p[1:len(self.t) + 1], kind=self.kind)
         t = self.get_t(p, dt, cutoff, maxtmax)
         s = p[0] * f(t)
diff --git a/pastas/solver.py b/pastas/solver.py
index fa0331ac8..59093420d 100644
--- a/pastas/solver.py
+++ b/pastas/solver.py
@@ -45,7 +45,12 @@ class BaseSolver:
 
     """
 
-    def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[int] = None, obj_func: Optional[Function] = None, **kwargs) -> None:
+    def __init__(self,
+                 ml: Model,
+                 pcov: Optional[DataFrame] = None,
+                 nfev: Optional[int] = None,
+                 obj_func: Optional[Function] = None,
+                 **kwargs) -> None:
         self.ml = ml
         self.pcov = pcov  # Covariances of the parameters
         if pcov is None:
@@ -56,8 +61,14 @@ def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[i
         self.obj_func = obj_func
         self.result = None  # Object returned by the optimization method
 
-    def misfit(self, p: ArrayLike, noise: bool, weights: Optional[Series] = None, callback: Optional[CallBack] = None,
-               returnseparate: bool = False) -> Union[ArrayLike, Tuple[ArrayLike, ArrayLike, ArrayLike]]:
+    def misfit(
+        self,
+        p: ArrayLike,
+        noise: bool,
+        weights: Optional[Series] = None,
+        callback: Optional[CallBack] = None,
+        returnseparate: bool = False
+    ) -> Union[ArrayLike, Tuple[ArrayLike, ArrayLike, ArrayLike]]:
         """This method is called by all solvers to obtain a series that are
         minimized in the optimization process. It handles the application of
         the weights, a noisemodel and other optimization options.
@@ -105,7 +116,11 @@ def misfit(self, p: ArrayLike, noise: bool, weights: Optional[Series] = None, ca
 
         return rv.values
 
-    def prediction_interval(self, n: int = 1000, alpha: float = 0.05, max_iter: int = 10, **kwargs) -> DataFrame:
+    def prediction_interval(self,
+                            n: int = 1000,
+                            alpha: float = 0.05,
+                            max_iter: int = 10,
+                            **kwargs) -> DataFrame:
         """Method to calculate the prediction interval for the simulation.
 
         Returns
@@ -131,7 +146,11 @@ def prediction_interval(self, n: int = 1000, alpha: float = 0.05, max_iter: int
         rv = data.quantile(q, axis=1).transpose()
         return rv
 
-    def ci_simulation(self, n: int = 1000, alpha: float = 0.05, max_iter: int = 10, **kwargs) -> DataFrame:
+    def ci_simulation(self,
+                      n: int = 1000,
+                      alpha: float = 0.05,
+                      max_iter: int = 10,
+                      **kwargs) -> DataFrame:
         """Method to calculate the confidence interval for the simulation.
 
         Returns
@@ -152,7 +171,12 @@ def ci_simulation(self, n: int = 1000, alpha: float = 0.05, max_iter: int = 10,
                                              alpha=alpha, max_iter=max_iter,
                                              **kwargs)
 
-    def ci_block_response(self, name: str, n: int = 1000, alpha: float = 0.05, max_iter: int = 10, **kwargs) -> DataFrame:
+    def ci_block_response(self,
+                          name: str,
+                          n: int = 1000,
+                          alpha: float = 0.05,
+                          max_iter: int = 10,
+                          **kwargs) -> DataFrame:
         """Method to calculate the confidence interval for the block response.
 
         Returns
@@ -175,7 +199,12 @@ def ci_block_response(self, name: str, n: int = 1000, alpha: float = 0.05, max_i
                                              max_iter=max_iter, dt=dt,
                                              **kwargs)
 
-    def ci_step_response(self, name: str, n: int = 1000, alpha: float = 0.05, max_iter: int = 10, **kwargs) -> DataFrame:
+    def ci_step_response(self,
+                         name: str,
+                         n: int = 1000,
+                         alpha: float = 0.05,
+                         max_iter: int = 10,
+                         **kwargs) -> DataFrame:
         """Method to calculate the confidence interval for the step response.
 
         Returns
@@ -198,7 +227,12 @@ def ci_step_response(self, name: str, n: int = 1000, alpha: float = 0.05, max_it
                                              max_iter=max_iter, dt=dt,
                                              **kwargs)
 
-    def ci_contribution(self, name: str, n: int = 1000, alpha: float = 0.05, max_iter: int = 10, **kwargs) -> DataFrame:
+    def ci_contribution(self,
+                        name: str,
+                        n: int = 1000,
+                        alpha: float = 0.05,
+                        max_iter: int = 10,
+                        **kwargs) -> DataFrame:
         """Method to calculate the confidence interval for the contribution.
 
         Returns
@@ -220,7 +254,10 @@ def ci_contribution(self, name: str, n: int = 1000, alpha: float = 0.05, max_ite
                                              max_iter=max_iter,
                                              **kwargs)
 
-    def get_parameter_sample(self, name: Optional[str] = None, n: int = None, max_iter: int = 10) -> ArrayLike:
+    def get_parameter_sample(self,
+                             name: Optional[str] = None,
+                             n: int = None,
+                             max_iter: int = 10) -> ArrayLike:
         """Method to obtain a parameter sets for monte carlo analyses.
 
         Parameters
@@ -279,7 +316,11 @@ def get_parameter_sample(self, name: Optional[str] = None, n: int = None, max_it
                            f"{samples.shape[0]}/{n}. Increase 'max_iter'.")
         return samples[:n, :]
 
-    def _get_realizations(self, func: Function, n: Optional[int] = None, name: Optional[str] = None, max_iter: int = 10,
+    def _get_realizations(self,
+                          func: Function,
+                          n: Optional[int] = None,
+                          name: Optional[str] = None,
+                          max_iter: int = 10,
                           **kwargs) -> DataFrame:
         """Internal method to obtain n number of parameter realizations."""
         if name:
@@ -294,7 +335,13 @@ def _get_realizations(self, func: Function, n: Optional[int] = None, name: Optio
 
         return DataFrame.from_dict(data, orient="columns", dtype=float)
 
-    def _get_confidence_interval(self,  func: Function, n: Optional[int] = None, name: Optional[str] = None, max_iter: int = 10, alpha: float = 0.05, **kwargs) -> DataFrame:
+    def _get_confidence_interval(self,
+                                 func: Function,
+                                 n: Optional[int] = None,
+                                 name: Optional[str] = None,
+                                 max_iter: int = 10,
+                                 alpha: float = 0.05,
+                                 **kwargs) -> DataFrame:
         """Internal method to obtain a confidence interval."""
         q = [alpha / 2, 1 - alpha / 2]
         data = self._get_realizations(func=func, n=n, name=name,
@@ -362,7 +409,11 @@ def to_dict(self) -> dict:
 class LeastSquares(BaseSolver):
     _name = "LeastSquares"
 
-    def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[int] = None, **kwargs) -> None:
+    def __init__(self,
+                 ml: Model,
+                 pcov: Optional[DataFrame] = None,
+                 nfev: Optional[int] = None,
+                 **kwargs) -> None:
         """Solver based on Scipy's least_squares method [scipy_ref]_.
 
         Notes
@@ -384,7 +435,11 @@ def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[i
         """
         BaseSolver.__init__(self, ml=ml, pcov=pcov, nfev=nfev, **kwargs)
 
-    def solve(self, noise: bool = True, weights: Optional[Series] = None, callback: Optional[CallBack] = None, **kwargs) -> Tuple[bool, ArrayLike, ArrayLike]:
+    def solve(self,
+              noise: bool = True,
+              weights: Optional[Series] = None,
+              callback: Optional[CallBack] = None,
+              **kwargs) -> Tuple[bool, ArrayLike, ArrayLike]:
         self.vary = self.ml.parameters.vary.values.astype(bool)
         self.initial = self.ml.parameters.initial.values.copy()
         parameters = self.ml.parameters.loc[self.vary]
@@ -413,13 +468,17 @@ def solve(self, noise: bool = True, weights: Optional[Series] = None, callback:
 
         return success, optimal, stderr
 
-    def objfunction(self, p: ArrayLike, noise: bool, weights: Series, callback: CallBack) -> ArrayLike:
+    def objfunction(self, p: ArrayLike, noise: bool, weights: Series,
+                    callback: CallBack) -> ArrayLike:
         par = self.initial
         par[self.vary] = p
         return self.misfit(p=par, noise=noise, weights=weights,
                            callback=callback)
 
-    def _get_covariances(self, jacobian: ArrayLike, cost: float, absolute_sigma: bool = False) -> ArrayLike:
+    def _get_covariances(self,
+                         jacobian: ArrayLike,
+                         cost: float,
+                         absolute_sigma: bool = False) -> ArrayLike:
         """Internal method to get the covariance matrix from the jacobian.
 
         Parameters
@@ -473,7 +532,11 @@ def _get_covariances(self, jacobian: ArrayLike, cost: float, absolute_sigma: boo
 class LmfitSolve(BaseSolver):
     _name = "LmfitSolve"
 
-    def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[int] = None, **kwargs) -> None:
+    def __init__(self,
+                 ml: Model,
+                 pcov: Optional[DataFrame] = None,
+                 nfev: Optional[int] = None,
+                 **kwargs) -> None:
         """Solving the model using the LmFit solver [LM]_.
 
          This is basically a wrapper around the scipy solvers, adding some
@@ -492,7 +555,11 @@ def __init__(self, ml: Model, pcov: Optional[DataFrame] = None, nfev: Optional[i
             raise ImportError(msg)
         BaseSolver.__init__(self, ml=ml, pcov=pcov, nfev=nfev, **kwargs)
 
-    def solve(self, noise: bool = True, weights: Optional[Series] = None, callback: Optional[CallBack] = None, method: Optional[str] = "leastsq",
+    def solve(self,
+              noise: bool = True,
+              weights: Optional[Series] = None,
+              callback: Optional[CallBack] = None,
+              method: Optional[str] = "leastsq",
               **kwargs) -> Tuple[bool, ArrayLike, ArrayLike]:
 
         # Deal with the parameters
@@ -536,6 +603,7 @@ def solve(self, noise: bool = True, weights: Optional[Series] = None, callback:
 
         return success, optimal[:idx], stderr[:idx]
 
-    def objfunction(self, parameters: DataFrame, noise: bool, weights: Series, callback: CallBack) -> ArrayLike:
+    def objfunction(self, parameters: DataFrame, noise: bool, weights: Series,
+                    callback: CallBack) -> ArrayLike:
         p = np.array([p.value for p in parameters.values()])
         return self.misfit(p=p, noise=noise, weights=weights, callback=callback)
diff --git a/pastas/stats/core.py b/pastas/stats/core.py
index 2f9e8f7b5..f623de8c2 100644
--- a/pastas/stats/core.py
+++ b/pastas/stats/core.py
@@ -18,8 +18,14 @@
 from pastas.typing import ArrayLike
 
 
-def acf(x: Series, lags: ArrayLike = 365, bin_method: str = 'rectangle', bin_width: float = 0.5, max_gap: float = inf,
-        min_obs: int = 20, full_output: bool = False, alpha: float = 0.05) -> Union[Series, DataFrame]:
+def acf(x: Series,
+        lags: ArrayLike = 365,
+        bin_method: str = 'rectangle',
+        bin_width: float = 0.5,
+        max_gap: float = inf,
+        min_obs: int = 20,
+        full_output: bool = False,
+        alpha: float = 0.05) -> Union[Series, DataFrame]:
     """Calculate the autocorrelation function for irregular time steps.
 
     Parameters
@@ -89,8 +95,15 @@ def acf(x: Series, lags: ArrayLike = 365, bin_method: str = 'rectangle', bin_wid
         return c
 
 
-def ccf(x: Series, y: Series, lags: ArrayLike = 365, bin_method: str = 'rectangle', bin_width: float = 0.5,
-        max_gap: float = inf, min_obs: int = 20, full_output: bool = False, alpha: float = 0.05) -> Union[Series, DataFrame]:
+def ccf(x: Series,
+        y: Series,
+        lags: ArrayLike = 365,
+        bin_method: str = 'rectangle',
+        bin_width: float = 0.5,
+        max_gap: float = inf,
+        min_obs: int = 20,
+        full_output: bool = False,
+        alpha: float = 0.05) -> Union[Series, DataFrame]:
     """Method to compute the cross-correlation for irregular time series.
 
     Parameters
@@ -194,7 +207,13 @@ def _preprocess(x: Series, max_gap: int) -> Tuple[Series, float, float]:
 
 
 @njit
-def _compute_ccf_rectangle(lags: ArrayLike, t_x: ArrayLike, x: ArrayLike, t_y: ArrayLike, y: ArrayLike, bin_width: float = 0.5) -> Tuple[ArrayLike, ArrayLike]:
+def _compute_ccf_rectangle(
+        lags: ArrayLike,
+        t_x: ArrayLike,
+        x: ArrayLike,
+        t_y: ArrayLike,
+        y: ArrayLike,
+        bin_width: float = 0.5) -> Tuple[ArrayLike, ArrayLike]:
     """Internal numba-optimized method to compute the ccf."""
     c = empty_like(lags)
     b = empty_like(lags)
@@ -220,7 +239,13 @@ def _compute_ccf_rectangle(lags: ArrayLike, t_x: ArrayLike, x: ArrayLike, t_y: A
 
 
 @njit
-def _compute_ccf_gaussian(lags: ArrayLike, t_x: ArrayLike, x: ArrayLike, t_y: ArrayLike, y: ArrayLike, bin_width: float = 0.5) -> Tuple[ArrayLike, ArrayLike]:
+def _compute_ccf_gaussian(
+        lags: ArrayLike,
+        t_x: ArrayLike,
+        x: ArrayLike,
+        t_y: ArrayLike,
+        y: ArrayLike,
+        bin_width: float = 0.5) -> Tuple[ArrayLike, ArrayLike]:
     """Internal numba-optimized method to compute the ccf."""
     c = empty_like(lags)
     b = empty_like(lags)
@@ -249,7 +274,8 @@ def _compute_ccf_gaussian(lags: ArrayLike, t_x: ArrayLike, x: ArrayLike, t_y: Ar
     return c, b
 
 
-def _compute_ccf_regular(lags: ArrayLike, x: ArrayLike, y: ArrayLike) -> Tuple[ArrayLike, ArrayLike]:
+def _compute_ccf_regular(lags: ArrayLike, x: ArrayLike,
+                         y: ArrayLike) -> Tuple[ArrayLike, ArrayLike]:
     c = empty_like(lags)
     for i, lag in enumerate(lags):
         c[i] = corrcoef(x[:-int(lag)], y[int(lag):])[0, 1]
diff --git a/pastas/stats/dutch.py b/pastas/stats/dutch.py
index 532c6bb9a..b9ab9dfde 100644
--- a/pastas/stats/dutch.py
+++ b/pastas/stats/dutch.py
@@ -13,9 +13,11 @@
 from pastas.typing import Function, TimestampType
 
 
-def q_ghg(
-        series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, q: float = 0.94,
-        by_year: bool = True) -> Series:
+def q_ghg(series: Series,
+          tmin: Optional[TimestampType] = None,
+          tmax: Optional[TimestampType] = None,
+          q: float = 0.94,
+          by_year: bool = True) -> Series:
     """Gemiddeld Hoogste Grondwaterstand (GHG) also called MHGL (Mean High
     Groundwater Level).
 
@@ -39,9 +41,11 @@ def q_ghg(
     return _q_gxg(series, q, tmin=tmin, tmax=tmax, by_year=by_year)
 
 
-def q_glg(
-        series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, q: float = 0.06,
-        by_year: bool = True) -> Series:
+def q_glg(series: Series,
+          tmin: Optional[TimestampType] = None,
+          tmax: Optional[TimestampType] = None,
+          q: float = 0.06,
+          by_year: bool = True) -> Series:
     """Gemiddeld Laagste Grondwaterstand (GLG) also called MLGL (Mean Low
     Groundwater Level).
 
@@ -65,7 +69,10 @@ def q_glg(
     return _q_gxg(series, q, tmin=tmin, tmax=tmax, by_year=by_year)
 
 
-def q_gvg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, by_year: bool = True) -> Series:
+def q_gvg(series: Series,
+          tmin: Optional[TimestampType] = None,
+          tmax: Optional[TimestampType] = None,
+          by_year: bool = True) -> Series:
     """Gemiddeld Voorjaarsgrondwaterstand (GVG) also called MSGL (Mean Spring
     Groundwater Level).
 
@@ -104,8 +111,14 @@ def q_gvg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[T
         return nan
 
 
-def ghg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-        fill_method: str = 'nearest', limit: int = 0, output: str = 'mean', min_n_meas: int = 16, min_n_years: int = 8,
+def ghg(series: Series,
+        tmin: Optional[TimestampType] = None,
+        tmax: Optional[TimestampType] = None,
+        fill_method: str = 'nearest',
+        limit: int = 0,
+        output: str = 'mean',
+        min_n_meas: int = 16,
+        min_n_years: int = 8,
         year_offset: str = 'a-mar') -> Union[Series, float]:
     """Calculate the 'Gemiddelde Hoogste Grondwaterstand' (Average High
     Groundwater Level)
@@ -176,8 +189,14 @@ def mean_high(s, min_n_meas):
                 year_offset=year_offset)
 
 
-def glg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-        fill_method: str = 'nearest', limit: int = 0, output: str = 'mean', min_n_meas: int = 16, min_n_years: int = 8,
+def glg(series: Series,
+        tmin: Optional[TimestampType] = None,
+        tmax: Optional[TimestampType] = None,
+        fill_method: str = 'nearest',
+        limit: int = 0,
+        output: str = 'mean',
+        min_n_meas: int = 16,
+        min_n_years: int = 8,
         year_offset: str = 'a-mar') -> Union[Series, float]:
     """Calculate the 'Gemiddelde Laagste Grondwaterstand' (Average Low
     Groundwater Level).
@@ -249,8 +268,15 @@ def mean_low(s, min_n_meas):
                 year_offset=year_offset)
 
 
-def gvg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, fill_method: str = 'linear', limit: int = 8,
-        output: str = 'mean', min_n_meas: int = 2, min_n_years: int = 8, year_offset: str = 'a') -> Union[Series, float]:
+def gvg(series: Series,
+        tmin: Optional[TimestampType] = None,
+        tmax: Optional[TimestampType] = None,
+        fill_method: str = 'linear',
+        limit: int = 8,
+        output: str = 'mean',
+        min_n_meas: int = 2,
+        min_n_years: int = 8,
+        year_offset: str = 'a') -> Union[Series, float]:
     """Calculate the 'Gemiddelde Voorjaars Grondwaterstand' (Average Spring
     Groundwater Level).
 
@@ -308,8 +334,15 @@ def _mean_spring(s, min_n_meas):
                 year_offset=year_offset)
 
 
-def gg(series: Series, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, fill_method: str = 'nearest', limit: int = 0,
-        output: str = 'mean', min_n_meas: int = 16, min_n_years: int = 8, year_offset: str = 'a-mar') -> Union[Series, float]:
+def gg(series: Series,
+       tmin: Optional[TimestampType] = None,
+       tmax: Optional[TimestampType] = None,
+       fill_method: str = 'nearest',
+       limit: int = 0,
+       output: str = 'mean',
+       min_n_meas: int = 16,
+       min_n_years: int = 8,
+       year_offset: str = 'a-mar') -> Union[Series, float]:
     """Calculate the 'Gemiddelde Grondwaterstand' (Average Groundwater Level)
 
     Parameters
@@ -406,9 +439,9 @@ def isinspring(x): return (((x.month == 3) and (x.day >= 14)) or
 
 
 def _gxg(series: Series, year_agg: Function, tmin: Optional[TimestampType],
-         tmax: Optional[TimestampType],
-         fill_method: str, limit: Union[int, None],
-         output: str, min_n_meas: int, min_n_years: int, year_offset: str) -> Union[Series, float]:
+         tmax: Optional[TimestampType], fill_method: str,
+         limit: Union[int, None], output: str, min_n_meas: int,
+         min_n_years: int, year_offset: str) -> Union[Series, float]:
     """Internal method for classic GXG statistics. Resampling the series to
     every 14th and 28th of the month. Taking the mean of aggregated values per
     year.
@@ -550,9 +583,11 @@ def _gxg(series: Series, year_agg: Function, tmin: Optional[TimestampType],
         raise (ValueError(msg))
 
 
-def _q_gxg(
-        series: Series, q: float, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-        by_year: bool = True) -> Series:
+def _q_gxg(series: Series,
+           q: float,
+           tmin: Optional[TimestampType] = None,
+           tmax: Optional[TimestampType] = None,
+           by_year: bool = True) -> Series:
     """Dutch groundwater statistics GHG and GLG approximated by taking
     quantiles of the timeseries values per year and taking the mean of the
     quantiles.
diff --git a/pastas/stats/metrics.py b/pastas/stats/metrics.py
index bd291c0af..4bc99308e 100644
--- a/pastas/stats/metrics.py
+++ b/pastas/stats/metrics.py
@@ -27,7 +27,12 @@
 
 # Absolute Error Metrics
 
-def mae(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", weighted: bool = False,
+
+def mae(obs: Optional[Series] = None,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        weighted: bool = False,
         max_gap: int = 30) -> float:
     """Compute the (weighted) Mean Absolute Error (MAE).
 
@@ -75,7 +80,11 @@ def mae(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optiona
     return (w * abs(res.to_numpy())).sum()
 
 
-def rmse(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", weighted: bool = False,
+def rmse(obs: Optional[Series] = None,
+         sim: Optional[Series] = None,
+         res: Optional[Series] = None,
+         missing: str = "drop",
+         weighted: bool = False,
          max_gap: int = 30) -> float:
     """Compute the (weighted) Root Mean Squared Error (RMSE).
 
@@ -122,7 +131,10 @@ def rmse(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Option
     return sqrt((w * res.to_numpy() ** 2).sum())
 
 
-def sse(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop") -> float:
+def sse(obs: Optional[Series] = None,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop") -> float:
     """Compute the Sum of the Squared Errors (SSE).
 
     Parameters
@@ -162,7 +174,11 @@ def sse(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optiona
 
 # Percentage Error Metrics
 
-def pearsonr(obs: Series, sim: Series, missing: str = "drop", weighted: bool = False,
+
+def pearsonr(obs: Series,
+             sim: Series,
+             missing: str = "drop",
+             weighted: bool = False,
              max_gap: int = 30) -> float:
     """Compute the (weighted) Pearson correlation (r).
 
@@ -215,7 +231,11 @@ def pearsonr(obs: Series, sim: Series, missing: str = "drop", weighted: bool = F
     return r
 
 
-def evp(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", weighted: bool = False,
+def evp(obs: Series,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        weighted: bool = False,
         max_gap: int = 30) -> float:
     """Compute the (weighted) Explained Variance Percentage (EVP).
 
@@ -270,7 +290,11 @@ def evp(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None,
                          var(obs, weighted=weighted, max_gap=max_gap))) * 100
 
 
-def nse(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", weighted: bool = False,
+def nse(obs: Series,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        weighted: bool = False,
         max_gap: int = 30) -> float:
     """Compute the (weighted) Nash-Sutcliffe Efficiency (NSE).
 
@@ -319,8 +343,13 @@ def nse(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None,
         (w * (obs.to_numpy() - mu) ** 2).sum()
 
 
-def rsq(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", weighted: bool = False,
-        max_gap: int = 30,  nparam: Optional[int] = None) -> float:
+def rsq(obs: Series,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        weighted: bool = False,
+        max_gap: int = 30,
+        nparam: Optional[int] = None) -> float:
     """Compute R-squared, possibly adjusted for the number of free parameters.
 
     Parameters
@@ -379,7 +408,11 @@ def rsq(obs: Series, sim: Optional[Series] = None, res: Optional[Series] = None,
         return 1.0 - rss / tss
 
 
-def bic(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", nparam: int = 1) -> float:
+def bic(obs: Optional[Series] = None,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        nparam: int = 1) -> float:
     """Compute the Bayesian Information Criterium (BIC).
 
     Parameters
@@ -423,7 +456,11 @@ def bic(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optiona
             nparam * log(res.index.size))
 
 
-def aic(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optional[Series] = None, missing: str = "drop", nparam: int = 1) -> float:
+def aic(obs: Optional[Series] = None,
+        sim: Optional[Series] = None,
+        res: Optional[Series] = None,
+        missing: str = "drop",
+        nparam: int = 1) -> float:
     """Compute the Akaike Information Criterium (AIC).
 
     Parameters
@@ -468,8 +505,10 @@ def aic(obs: Optional[Series] = None, sim: Optional[Series] = None, res: Optiona
 
 
 # Forecast Error Metrics
-
-def kge_2012(obs: Series, sim: Series, missing: str = "drop", weighted: bool = False,
+def kge_2012(obs: Series,
+             sim: Series,
+             missing: str = "drop",
+             weighted: bool = False,
              max_gap: int = 30) -> float:
     """Compute the (weighted) Kling-Gupta Efficiency (KGE).
 
diff --git a/pastas/stats/signatures.py b/pastas/stats/signatures.py
index 133367c60..fdb81bd6c 100644
--- a/pastas/stats/signatures.py
+++ b/pastas/stats/signatures.py
@@ -263,7 +263,10 @@ def interannual_variation(series: Series, normalize: bool = True) -> float:
     return (hw.max() - hw.min()) + (hl.max() - hl.min()) / 2
 
 
-def colwell_components(series: Series, bins: int = 11, freq: str = "M", method: str = "mean",
+def colwell_components(series: Series,
+                       bins: int = 11,
+                       freq: str = "M",
+                       method: str = "mean",
                        normalize: bool = True) -> Tuple[float, float, float]:
     """Colwell predictability, constant, and contingency [colwell_1974]_.
 
@@ -333,7 +336,10 @@ def colwell_components(series: Series, bins: int = 11, freq: str = "M", method:
     return p, c, m
 
 
-def colwell_constancy(series: Series, bins: int = 11, freq: str = "M", method: str = "mean",
+def colwell_constancy(series: Series,
+                      bins: int = 11,
+                      freq: str = "M",
+                      method: str = "mean",
                       normalize: bool = True) -> Tuple[float, float, float]:
     """Colwells constancy index after [colwell_1974]_.
 
@@ -366,12 +372,17 @@ def colwell_constancy(series: Series, bins: int = 11, freq: str = "M", method: s
        Contingency of periodic phenomena. Ecology, 55(5), 1148–1153.
 
     """
-    return \
-        colwell_components(series=series, bins=bins, freq=freq, method=method,
-                           normalize=normalize)[1]
+    return colwell_components(series=series,
+                              bins=bins,
+                              freq=freq,
+                              method=method,
+                              normalize=normalize)[1]
 
 
-def colwell_contingency(series: Series, bins: int = 11, freq: str = "M", method: str = "mean",
+def colwell_contingency(series: Series,
+                        bins: int = 11,
+                        freq: str = "M",
+                        method: str = "mean",
                         normalize: bool = True) -> Tuple[float, float, float]:
     """Colwell contingency [colwell_1974]_
 
@@ -407,9 +418,11 @@ def colwell_contingency(series: Series, bins: int = 11, freq: str = "M", method:
        Contingency of periodic phenomena. Ecology, 55(5), 1148–1153.
 
     """
-    return \
-        colwell_components(series=series, bins=bins, freq=freq, method=method,
-                           normalize=normalize)[2]
+    return colwell_components(series=series,
+                              bins=bins,
+                              freq=freq,
+                              method=method,
+                              normalize=normalize)[2]
 
 
 def low_pulse_count(series: Series, quantile: float = 0.2) -> int:
@@ -855,7 +868,9 @@ def bimodality_coefficient(series: Series, normalize: bool = False) -> float:
            (kurt + ((3 * ((n - 1) ** 2)) / ((n - 2) * (n - 3))))
 
 
-def recession_constant(series: Series, bins: int = 20, normalize: bool = False) -> float:
+def recession_constant(series: Series,
+                       bins: int = 20,
+                       normalize: bool = False) -> float:
     """Recession constant after [kirchner_2009]_.
 
     Parameters
@@ -900,7 +915,9 @@ def recession_constant(series: Series, bins: int = 20, normalize: bool = False)
     return fit.slope
 
 
-def recovery_constant(series: Series, bins: int = 20, normalize: bool = False) -> float:
+def recovery_constant(series: Series,
+                      bins: int = 20,
+                      normalize: bool = False) -> float:
     """Recovery constant after [kirchner_2009]_.
 
     Parameters
@@ -946,7 +963,10 @@ def recovery_constant(series: Series, bins: int = 20, normalize: bool = False) -
     return fit.slope
 
 
-def duration_curve_slope(series: Series, l: float = 0.1, u: float = 0.9, normalize: bool = True) -> float:
+def duration_curve_slope(series: Series,
+                         l: float = 0.1,
+                         u: float = 0.9,
+                         normalize: bool = True) -> float:
     """Slope of the duration curve between percentile l and u.
 
     Parameters
@@ -985,7 +1005,10 @@ def duration_curve_slope(series: Series, l: float = 0.1, u: float = 0.9, normali
     return linregress(s.index, s.values).slope
 
 
-def duration_curve_range(series: Series, l: float = 0.1, u: float = 0.9, normalize: bool = True) -> float:
+def duration_curve_range(series: Series,
+                         l: float = 0.1,
+                         u: float = 0.9,
+                         normalize: bool = True) -> float:
     """Range of the duration curve between the percentile l and u.
 
     Parameters
diff --git a/pastas/stats/tests.py b/pastas/stats/tests.py
index 05a8660c5..32e7c980a 100644
--- a/pastas/stats/tests.py
+++ b/pastas/stats/tests.py
@@ -90,7 +90,10 @@ def durbin_watson(series: Series) -> float:
     return dw_stat, p
 
 
-def ljung_box(series: Series, lags: int = 15, nparam: int = 0, full_output: bool = False) -> Tuple[float, float]:
+def ljung_box(series: Series,
+              lags: int = 15,
+              nparam: int = 0,
+              full_output: bool = False) -> Tuple[float, float]:
     """Ljung-box test for autocorrelation.
 
     Parameters
@@ -278,8 +281,12 @@ def runs_test(series: Series, cutoff: str = "median") -> Tuple[float, float]:
     return z_stat, pval
 
 
-def stoffer_toloi(series: Series, lags: int = 15, nparam: int = 0, freq: str = "D",
-                  snap_to_equidistant_timestamps: bool = False) -> Tuple[float, float]:
+def stoffer_toloi(
+        series: Series,
+        lags: int = 15,
+        nparam: int = 0,
+        freq: str = "D",
+        snap_to_equidistant_timestamps: bool = False) -> Tuple[float, float]:
     """Adapted Ljung-Box test to deal with missing data [stoffer_1992]_.
 
     Parameters
@@ -415,7 +422,11 @@ def stoffer_toloi(series: Series, lags: int = 15, nparam: int = 0, freq: str = "
     return qm, pval
 
 
-def diagnostics(series: Series, alpha: float = 0.05, nparam: int = 0, lags: int = 15, stats: tuple = (),
+def diagnostics(series: Series,
+                alpha: float = 0.05,
+                nparam: int = 0,
+                lags: int = 15,
+                stats: tuple = (),
                 float_fmt: str = "{0:.2f}") -> DataFrame:
     """Methods to compute various diagnostics checks for a time series.
 
diff --git a/pastas/stressmodels.py b/pastas/stressmodels.py
index a1dcfa2bb..11466de2c 100644
--- a/pastas/stressmodels.py
+++ b/pastas/stressmodels.py
@@ -51,9 +51,14 @@ class StressModelBase:
     """
     _name = "StressModelBase"
 
-    def __init__(
-            self, name: str, tmin: TimestampType, tmax: TimestampType, rfunc: Optional[RFunc] = None, up: bool = True,
-            meanstress: float = 1.0, cutoff: float = 0.999) -> None:
+    def __init__(self,
+                 name: str,
+                 tmin: TimestampType,
+                 tmax: TimestampType,
+                 rfunc: Optional[RFunc] = None,
+                 up: bool = True,
+                 meanstress: float = 1.0,
+                 cutoff: float = 0.999) -> None:
         self.name = validate_name(name)
         self.tmin = tmin
         self.tmax = tmax
@@ -162,9 +167,13 @@ def dump_stress(self, series: bool = True) -> dict:
 
         return data
 
-    def get_stress(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, istress: Optional[int] = None, **kwargs) -> DataFrame:
+    def get_stress(self,
+                   p: Optional[ArrayLike] = None,
+                   tmin: Optional[TimestampType] = None,
+                   tmax: Optional[TimestampType] = None,
+                   freq: Optional[str] = None,
+                   istress: Optional[int] = None,
+                   **kwargs) -> DataFrame:
         """Returns the stress or stresses of the time series object as a pandas
         DataFrame.
 
@@ -208,7 +217,8 @@ def get_nsplit(self) -> int:
         else:
             return len(self.stress)
 
-    def _get_block(self, p: ArrayLike, dt: float, tmin: TimestampType, tmax: TimestampType) -> ArrayLike:
+    def _get_block(self, p: ArrayLike, dt: float, tmin: TimestampType,
+                   tmax: TimestampType) -> ArrayLike:
         """Internal method to get the block-response function."""
         if tmin is not None and tmax is not None:
             day = Timedelta(1, 'D')
@@ -286,9 +296,12 @@ def set_init_parameters(self) -> None:
         """Set the initial parameters (back) to their default values."""
         self.parameters = self.rfunc.get_init_parameters(self.name)
 
-    def simulate(
-            self, p: ArrayLike, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0) -> Series:
+    def simulate(self,
+                 p: ArrayLike,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0) -> Series:
         """Simulates the head contribution.
 
         Parameters
@@ -364,7 +377,12 @@ class StepModel(StressModelBase):
     """
     _name = "StepModel"
 
-    def __init__(self, tstart: TimestampType, name: str, rfunc: Optional[RFunc] = None, up: bool = True, cutoff: float = 0.999) -> None:
+    def __init__(self,
+                 tstart: TimestampType,
+                 name: str,
+                 rfunc: Optional[RFunc] = None,
+                 up: bool = True,
+                 cutoff: float = 0.999) -> None:
         if rfunc is None:
             rfunc = One()
         StressModelBase.__init__(self, name=name, tmin=Timestamp.min,
@@ -382,9 +400,12 @@ def set_init_parameters(self) -> None:
         self.parameters.loc[self.name + "_tstart"] = (tinit, tmin, tmax,
                                                       False, self.name)
 
-    def simulate(
-            self, p: ArrayLike, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0) -> Series:
+    def simulate(self,
+                 p: ArrayLike,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0) -> Series:
         tstart = Timestamp.fromordinal(int(p[-1]))
         tindex = date_range(tmin, tmax, freq=freq)
         h = Series(0, tindex, name=self.name)
@@ -429,8 +450,13 @@ class LinearTrend(StressModelBase):
     """
     _name = "LinearTrend"
 
-    def __init__(self, start: TimestampType, end: TimestampType, name: str = "trend") -> None:
-        StressModelBase.__init__(self, name=name, tmin=Timestamp.min,
+    def __init__(self,
+                 start: TimestampType,
+                 end: TimestampType,
+                 name: str = "trend") -> None:
+        StressModelBase.__init__(self,
+                                 name=name,
+                                 tmin=Timestamp.min,
                                  tmax=Timestamp.max)
         self.start = start
         self.end = end
@@ -450,9 +476,12 @@ def set_init_parameters(self) -> None:
         self.parameters.loc[self.name + "_tend"] = (end, tmin, tmax,
                                                     False, self.name)
 
-    def simulate(
-            self, p: ArrayLike, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0) -> Series:
+    def simulate(self,
+                 p: ArrayLike,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0) -> Series:
         """Simulate the trend."""
         tindex = date_range(tmin, tmax, freq=freq)
 
@@ -550,10 +579,15 @@ class WellModel(StressModelBase):
     """
     _name = "WellModel"
 
-    def __init__(
-            self, stress: List[Series],
-            rfunc: RFunc, name: str, distances: ArrayLike, up: bool = False, cutoff: float = 0.999, settings: str = "well",
-            sort_wells: bool = True) -> None:
+    def __init__(self,
+                 stress: List[Series],
+                 rfunc: RFunc,
+                 name: str,
+                 distances: ArrayLike,
+                 up: bool = False,
+                 cutoff: float = 0.999,
+                 settings: str = "well",
+                 sort_wells: bool = True) -> None:
         if not (isinstance(rfunc, HantushWellModel) or
                 issubclass(rfunc, HantushWellModel)):
             raise NotImplementedError("WellModel only supports the rfunc "
@@ -612,9 +646,14 @@ def __init__(
     def set_init_parameters(self) -> None:
         self.parameters = self.rfunc.get_init_parameters(self.name)
 
-    def simulate(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0, istress: Optional[int] = None, **kwargs) -> Series:
+    def simulate(self,
+                 p: Optional[ArrayLike] = None,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0,
+                 istress: Optional[int] = None,
+                 **kwargs) -> Series:
         distances = self.get_distances(istress=istress)
         stress_df = self.get_stress(p=p, tmin=tmin, tmax=tmax, freq=freq,
                                     istress=istress, squeeze=False)
@@ -669,9 +708,14 @@ def handle_stress(stress, settings):
                          "dict or list.")
         return data
 
-    def get_stress(self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None,
-                   tmax: Optional[TimestampType] = None, freq: Optional[str] = None, istress: Optional[int] = None,
-                   squeeze: bool = True, **kwargs) -> DataFrame:
+    def get_stress(self,
+                   p: Optional[ArrayLike] = None,
+                   tmin: Optional[TimestampType] = None,
+                   tmax: Optional[TimestampType] = None,
+                   freq: Optional[str] = None,
+                   istress: Optional[int] = None,
+                   squeeze: bool = True,
+                   **kwargs) -> DataFrame:
         if tmin is None:
             tmin = self.tmin
         if tmax is None:
@@ -758,7 +802,10 @@ def to_dict(self, series: bool = True) -> dict:
         }
         return data
 
-    def variance_gain(self, model: Model, istress: Optional[int] = None, r: Optional[ArrayLike] = None) -> float:
+    def variance_gain(self,
+                      model: Model,
+                      istress: Optional[int] = None,
+                      r: Optional[ArrayLike] = None) -> float:
         """Calculate variance of the gain for WellModel.
 
         Variance of the gain is calculated based on propagation of uncertainty
@@ -880,12 +927,18 @@ class RechargeModel(StressModelBase):
     _name = "RechargeModel"
 
     def __init__(
-            self, prec: Series, evap: Series, rfunc: Optional[RFunc] = None, name: str = "recharge",
-            recharge: Optional[Recharge] = None, temp: Series = None, cutoff: float = 0.999,
-            settings: Tuple[Union[str, dict],
-                            Union[str, dict],
-                            Union[str, dict]] = ("prec", "evap", "evap"),
-            metadata: Optional[Tuple[dict, dict, dict]] = (None, None, None)) -> None:
+        self,
+        prec: Series,
+        evap: Series,
+        rfunc: Optional[RFunc] = None,
+        name: str = "recharge",
+        recharge: Optional[Recharge] = None,
+        temp: Series = None,
+        cutoff: float = 0.999,
+        settings: Tuple[Union[str, dict], Union[str, dict],
+                        Union[str, dict]] = ("prec", "evap", "evap"),
+        metadata: Optional[Tuple[dict, dict, dict]] = (None, None, None)
+    ) -> None:
 
         if rfunc is None:
             rfunc = Exponential()
@@ -985,9 +1038,14 @@ def update_stress(self, **kwargs) -> None:
         if "freq" in kwargs:
             self.freq = kwargs["freq"]
 
-    def simulate(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0, istress: Optional[int] = None, **kwargs) -> Series:
+    def simulate(self,
+                 p: Optional[ArrayLike] = None,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0,
+                 istress: Optional[int] = None,
+                 **kwargs) -> Series:
         """Method to simulate the contribution of recharge to the head.
 
         Parameters
@@ -1024,9 +1082,13 @@ def simulate(
         return Series(data=fftconvolve(stress, b, 'full')[:stress.size],
                       index=self.prec.series.index, name=name, fastpath=True)
 
-    def get_stress(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, istress: Optional[int] = None, **kwargs) -> Series:
+    def get_stress(self,
+                   p: Optional[ArrayLike] = None,
+                   tmin: Optional[TimestampType] = None,
+                   tmax: Optional[TimestampType] = None,
+                   freq: Optional[str] = None,
+                   istress: Optional[int] = None,
+                   **kwargs) -> Series:
         """Method to obtain the recharge stress calculated by the model.
 
         Parameters
@@ -1075,9 +1137,11 @@ def get_stress(
         else:
             return self.temp.series
 
-    def get_water_balance(
-            self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None) -> DataFrame:
+    def get_water_balance(self,
+                          p: Optional[ArrayLike] = None,
+                          tmin: Optional[TimestampType] = None,
+                          tmax: Optional[TimestampType] = None,
+                          freq: Optional[str] = None) -> DataFrame:
         """Method to obtain the water balance components.
 
         Parameters
@@ -1191,9 +1255,14 @@ class TarsoModel(RechargeModel):
     """
     _name = "TarsoModel"
 
-    def __init__(
-            self, prec: Series, evap: Series, oseries: Optional[Series] = None, dmin: Optional[float] = None,
-            dmax: Optional[float] = None, rfunc: Optional[RFunc] = Exponential, **kwargs) -> None:
+    def __init__(self,
+                 prec: Series,
+                 evap: Series,
+                 oseries: Optional[Series] = None,
+                 dmin: Optional[float] = None,
+                 dmax: Optional[float] = None,
+                 rfunc: Optional[RFunc] = Exponential,
+                 **kwargs) -> None:
         check_numba()
         if oseries is not None:
             if dmin is not None or dmax is not None:
@@ -1239,8 +1308,12 @@ def set_init_parameters(self) -> None:
         # combine all parameters
         self.parameters = concat([p0, p1, pr])
 
-    def simulate(self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None,
-                 tmax: Optional[TimestampType] = None, freq=None, dt: float = 1.0) -> Series:
+    def simulate(self,
+                 p: Optional[ArrayLike] = None,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq=None,
+                 dt: float = 1.0) -> Series:
         stress = self.get_stress(p=p, tmin=tmin, tmax=tmax, freq=freq)
         h = self.tarso(p[:-self.recharge.nparam], stress.values, dt)
         sim = Series(h, name=self.name, index=stress.index)
@@ -1352,8 +1425,16 @@ class ChangeModel(StressModelBase):
     """
     _name = "ChangeModel"
 
-    def __init__(self, stress: Series, rfunc1: RFunc, rfunc2: RFunc, name: str, tchange: str, up: bool = True,
-                 cutoff: float = 0.999, settings: Optional[Union[dict, str]] = None, metadata: Optional[dict] = None) -> None:
+    def __init__(self,
+                 stress: Series,
+                 rfunc1: RFunc,
+                 rfunc2: RFunc,
+                 name: str,
+                 tchange: str,
+                 up: bool = True,
+                 cutoff: float = 0.999,
+                 settings: Optional[Union[dict, str]] = None,
+                 metadata: Optional[dict] = None) -> None:
         if isinstance(stress, list):
             stress = stress[0]  # TODO Temporary fix Raoul, 2017-10-24
 
@@ -1400,9 +1481,12 @@ def set_init_parameters(self) -> None:
                                                        False, self.name)
         self.parameters.name = self.name
 
-    def simulate(
-            self, p: ArrayLike, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0) -> Series:
+    def simulate(self,
+                 p: ArrayLike,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0) -> Series:
         self.update_stress(tmin=tmin, tmax=tmax, freq=freq)
         rfunc1 = self.rfunc1.block(p[:self.rfunc1.nparam])
         rfunc2 = self.rfunc2.block(
@@ -1460,12 +1544,15 @@ class ReservoirModel(StressModelBase):
     """
     _name = "ReservoirModel"
 
-    def __init__(
-            self, stress: Series, reservoir: Reservoir, name: str, meanhead: float,
-            settings: Optional[Tuple[Union[str, dict],
-                                     Union[str, dict]]] = ("prec", "evap"),
-            metadata: Optional[Tuple[dict, dict]] = (None, None),
-            meanstress: Optional[float] = None) -> None:
+    def __init__(self,
+                 stress: Series,
+                 reservoir: Reservoir,
+                 name: str,
+                 meanhead: float,
+                 settings: Optional[Tuple[Union[str, dict],
+                                          Union[str, dict]]] = ("prec", "evap"),
+                 metadata: Optional[Tuple[dict, dict]] = (None, None),
+                 meanstress: Optional[float] = None) -> None:
         # Set resevoir object
         self.reservoir = reservoir(meanhead)
 
@@ -1504,9 +1591,13 @@ def set_init_parameters(self) -> None:
         """Set the initial parameters back to their default values."""
         self.parameters = self.reservoir.get_init_parameters(self.name)
 
-    def simulate(
-            self, p: ArrayLike, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None,
-            freq: Optional[str] = None, dt: float = 1.0, istress: Optional[int] = None) -> Series:
+    def simulate(self,
+                 p: ArrayLike,
+                 tmin: Optional[TimestampType] = None,
+                 tmax: Optional[TimestampType] = None,
+                 freq: Optional[str] = None,
+                 dt: float = 1.0,
+                 istress: Optional[int] = None) -> Series:
         """Simulates the head contribution.
 
         Parameters
@@ -1531,8 +1622,13 @@ def simulate(
                    index=stress[0].index, name=self.name, fastpath=True)
         return h
 
-    def get_stress(self, p: Optional[ArrayLike] = None, tmin: Optional[TimestampType] = None, tmax: Optional[TimestampType] = None, freq: Optional[str] = None,
-                   istress: int = 0, **kwargs) -> Tuple[Series, Series]:
+    def get_stress(self,
+                   p: Optional[ArrayLike] = None,
+                   tmin: Optional[TimestampType] = None,
+                   tmax: Optional[TimestampType] = None,
+                   freq: Optional[str] = None,
+                   istress: int = 0,
+                   **kwargs) -> Tuple[Series, Series]:
         if tmin is None:
             tmin = self.tmin
         if tmax is None:
diff --git a/pastas/timer.py b/pastas/timer.py
index e983802af..2a31eff83 100644
--- a/pastas/timer.py
+++ b/pastas/timer.py
@@ -40,7 +40,10 @@ class SolveTimer(tqdm):
     be updated quite as nicely.
     """
 
-    def __init__(self, max_time: Optional[float] = None, *args, **kwargs) -> None:
+    def __init__(self,
+                 max_time: Optional[float] = None,
+                 *args,
+                 **kwargs) -> None:
         """Initialize SolveTimer.
 
         Parameters
diff --git a/pastas/timeseries.py b/pastas/timeseries.py
index 3751a3634..2cbe5c311 100644
--- a/pastas/timeseries.py
+++ b/pastas/timeseries.py
@@ -60,8 +60,13 @@ class TimeSeries:
     """
     _predefined_settings = rcParams["timeseries"]
 
-    def __init__(self, series: Series, name: Optional[str] = None, settings: Optional[Union[str, dict]] = None, metadata: Optional[dict] = None,
-                 freq_original: str = None, **kwargs) -> None:
+    def __init__(self,
+                 series: Series,
+                 name: Optional[str] = None,
+                 settings: Optional[Union[str, dict]] = None,
+                 metadata: Optional[dict] = None,
+                 freq_original: str = None,
+                 **kwargs) -> None:
         if isinstance(series, TimeSeries):
             # Copy all the series
             self._series_original = series.series_original.copy()
diff --git a/pastas/transform.py b/pastas/transform.py
index 265473a74..5b4ccfb98 100644
--- a/pastas/transform.py
+++ b/pastas/transform.py
@@ -41,8 +41,12 @@ class ThresholdTransform:
     """
     _name = "ThresholdTransform"
 
-    def __init__(self, value: float = np.nan, vmin: float = np.nan, vmax: float = np.nan,
-                 name: str = 'ThresholdTransform', nparam: int = 2) -> None:
+    def __init__(self,
+                 value: float = np.nan,
+                 vmin: float = np.nan,
+                 vmax: float = np.nan,
+                 name: str = 'ThresholdTransform',
+                 nparam: int = 2) -> None:
         self.value = value
         self.vmin = vmin
         self.vmax = vmax
diff --git a/pastas/utils.py b/pastas/utils.py
index 5b25dd07b..d788ab7ea 100644
--- a/pastas/utils.py
+++ b/pastas/utils.py
@@ -301,7 +301,9 @@ def timestep_weighted_resample_fast(series0: Series, freq: str) -> Series:
     return series
 
 
-def get_equidistant_series(series: Series, freq: str, minimize_data_loss: bool = False) -> Series:
+def get_equidistant_series(series: Series,
+                           freq: str,
+                           minimize_data_loss: bool = False) -> Series:
     """Get equidistant timeseries using nearest reindexing.
 
     This method will shift observations to the nearest equidistant timestep to
@@ -474,7 +476,8 @@ def get_stress_tmin_tmax(ml: ModelType) -> Tuple[TimestampType, TimestampType]:
     return tmin, tmax
 
 
-def initialize_logger(logger: Optional[Any] = None, level: Optional[Any] = logging.INFO) -> None:
+def initialize_logger(logger: Optional[Any] = None,
+                      level: Optional[Any] = logging.INFO) -> None:
     """Internal method to create a logger instance to log program output.
 
     Parameters
@@ -492,7 +495,8 @@ def initialize_logger(logger: Optional[Any] = None, level: Optional[Any] = loggi
     # add_file_handlers(logger)
 
 
-def set_console_handler(logger: Optional[Any] = None, level: Optional[Any] = logging.INFO,
+def set_console_handler(logger: Optional[Any] = None,
+                        level: Optional[Any] = logging.INFO,
                         fmt: str = "%(levelname)s: %(message)s") -> None:
     """Method to add a console handler to the logger of Pastas.
 
@@ -549,11 +553,15 @@ def remove_console_handler(logger: Optional[Any] = None) -> None:
             logger.removeHandler(handler)
 
 
-def add_file_handlers(logger: Optional[Any] = None, filenames: Tuple[str] = ('info.log', 'errors.log'),
-                      levels: Tuple[Any] = (logging.INFO, logging.ERROR), maxBytes: int = 10485760,
-                      backupCount: int = 20, encoding: str = 'utf8',
-                      fmt: str = '%(asctime)s - %(name)s - %(levelname)s - %(message)s',
-                      datefmt: str = '%y-%m-%d %H:%M') -> None:
+def add_file_handlers(
+        logger: Optional[Any] = None,
+        filenames: Tuple[str] = ('info.log', 'errors.log'),
+        levels: Tuple[Any] = (logging.INFO, logging.ERROR),
+        maxBytes: int = 10485760,
+        backupCount: int = 20,
+        encoding: str = 'utf8',
+        fmt: str = '%(asctime)s - %(name)s - %(levelname)s - %(message)s',
+        datefmt: str = '%y-%m-%d %H:%M') -> None:
     """Method to add file handlers in the logger of Pastas.
 
     Parameters
@@ -620,7 +628,8 @@ def validate_name(name: str, raise_error: bool = False) -> str:
     name = str(name)
     for char in ilchar:
         if char in name:
-            msg = f"User-provided name '{name}' contains illegal character. Please remove {char} from name."
+            msg = f"User-provided name '{name}' contains illegal character."
+            msg += f"Please remove {char} from name."
             if raise_error:
                 raise Exception(msg)
             else: