docs: add comprehensive toolbox documentation

The documentation was sparse and consisted only of README with
installation instructions. This adds proper user guides, tutorials,
and example plots to help scientists get started with the library.

- Add quickstart guide with minimal examples
- Add batch reactor guide covering all Batch class functionality
- Add column transport guide covering all Column class functionality
- Add tutorial: organic matter degradation (Roden 2008)
- Add tutorial: sediment diagenesis with complex reaction network
- Add GitHub Pages workflow for automatic deployment
- Extract example plots from Jupyter notebooks
- Remove outdated publishing.md

Author: pycckuu

.github/workflows/docs.yml

@@ -0,0 +1,41 @@
+name: Deploy Documentation to GitHub Pages
+
+on:
+  push:
+    branches:
+      - master
+    paths:
+      - 'docs/**'
+      - '.github/workflows/docs.yml'
+  workflow_dispatch:
+
+permissions:
+  contents: read
+  pages: write
+  id-token: write
+
+concurrency:
+  group: "pages"
+  cancel-in-progress: false
+
+jobs:
+  deploy:
+    environment:
+      name: github-pages
+      url: ${{ steps.deployment.outputs.page_url }}
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v4
+
+      - name: Setup Pages
+        uses: actions/configure-pages@v4
+
+      - name: Upload artifact
+        uses: actions/upload-pages-artifact@v3
+        with:
+          path: 'docs'
+
+      - name: Deploy to GitHub Pages
+        id: deployment
+        uses: actions/deploy-pages@v4

docs/README.md

@@ -0,0 +1,52 @@
+# PorousMediaLab Documentation
+
+PorousMediaLab is a Python toolbox for batch and 1D reactive transport modeling in porous media, designed for scientists without extensive computational backgrounds.
+
+## Quick Links
+
+| Document | Description |
+|----------|-------------|
+| [Quickstart](quickstart.md) | Get started in 5 minutes |
+| [Batch Guide](batch-guide.md) | Complete guide to 0D batch reactor simulations |
+| [Column Guide](column-guide.md) | Complete guide to 1D reactive transport |
+
+## Tutorials
+
+Step-by-step tutorials with real-world examples:
+
+- [Organic Matter Degradation](tutorials/batch-roden.md) - Batch reactor modeling following Roden 2008
+- [Sediment Diagenesis](tutorials/sediment-column.md) - 1D column transport with complex reactions
+
+## Features
+
+- **Batch Module** - 0D closed-system reaction simulations
+- **Column Module** - 1D advection-diffusion-reaction equation solver
+- **Flexible Reactions** - Define rates using string expressions with numexpr
+- **pH Equilibrium** - Built-in acid-base chemistry calculations
+- **Henry Equilibrium** - Gas-liquid partitioning
+- **Calibration** - Parameter optimization against experimental data
+- **Visualization** - Built-in plotting methods
+
+## Installation
+
+```bash
+pip install porousmedialab
+```
+
+For development:
+
+```bash
+git clone https://github.com/biogeochemistry/PorousMediaLab
+cd PorousMediaLab
+poetry install
+```
+
+## Citation
+
+Igor Markelov (2020). Modelling Biogeochemical Cycles Across Scales: From Whole-Lake Phosphorus Dynamics to Microbial Reaction Systems. UWSpace. http://hdl.handle.net/10012/15513
+
+## Links
+
+- [GitHub Repository](https://github.com/biogeochemistry/PorousMediaLab)
+- [PyPI Package](https://pypi.org/project/porousmedialab/)
+- [Example Notebooks](https://github.com/biogeochemistry/PorousMediaLab/tree/master/examples)

docs/batch-guide.md

@@ -0,0 +1,212 @@
+# Batch Reactor Guide
+
+The `Batch` class simulates 0-dimensional (0D) batch reactor experiments - closed systems where only reactions occur without spatial transport.
+
+## Creating a Batch
+
+```python
+from porousmedialab.batch import Batch
+
+batch = Batch(tend=40, dt=1)
+```
+
+**Parameters:**
+- `tend` (float): End time of simulation (in your chosen time units)
+- `dt` (float): Timestep for integration
+
+## Adding Species
+
+```python
+batch.add_species(name='O2', init_conc=0.2)
+batch.add_species(name='CH4', init_conc=0.001)
+batch.add_species(name='CO2', init_conc=0)
+```
+
+**Parameters:**
+- `name` (str): Name of the chemical species (used in rate expressions)
+- `init_conc` (float): Initial concentration
+
+## Defining Constants
+
+Constants are parameters used in rate expressions:
+
+```python
+batch.constants['k1'] = 0.1          # Rate constant
+batch.constants['Km'] = 0.001        # Half-saturation constant
+batch.constants['T'] = 25            # Temperature
+```
+
+## Defining Rate Expressions
+
+Rate expressions are strings evaluated using [numexpr](https://github.com/pydata/numexpr):
+
+```python
+# First-order kinetics
+batch.rates['R1'] = 'k1 * A'
+
+# Michaelis-Menten kinetics
+batch.rates['R2'] = 'Vmax * S / (Km + S)'
+
+# Multiple terms
+batch.rates['R3'] = 'k1 * A * B / (Km + B)'
+
+# Inhibition term
+batch.rates['R4'] = 'k1 * A * Km_inh / (Km_inh + I)'
+```
+
+**Available operators:** `+`, `-`, `*`, `/`, `**` (power)
+
+**Available functions:** `exp`, `log`, `log10`, `sqrt`, `sin`, `cos`, `tan`, `abs`, `where`
+
+## Setting dcdt (Time Derivatives)
+
+Define how each species concentration changes over time:
+
+```python
+# Simple consumption and production
+batch.dcdt['A'] = '-R1'
+batch.dcdt['B'] = 'R1'
+
+# Multiple reactions affecting one species
+batch.dcdt['C'] = 'R1 - R2 + 0.5 * R3'
+
+# Stoichiometric coefficients
+batch.dcdt['O2'] = '-2 * R1 - 4 * R2'
+```
+
+## Running the Simulation
+
+```python
+batch.solve(verbose=True)
+```
+
+**Parameters:**
+- `verbose` (bool): If `True`, prints progress and time estimates
+
+## Accessing Results
+
+Results are stored in species dictionaries:
+
+```python
+# Get concentration time series (shape: [1, num_timesteps])
+conc = batch.A.concentration
+
+# Get final concentration
+final = batch.A.concentration[0, -1]
+
+# Get time array
+time = batch.time
+
+# Plot concentration vs time
+import matplotlib.pyplot as plt
+plt.plot(batch.time, batch.A.concentration[0])
+```
+
+## Plotting Methods
+
+Built-in plotting methods:
+
+```python
+# Plot single species
+batch.plot('A')
+
+# Plot all species
+batch.plot_profiles()
+
+# Plot reaction rates
+batch.plot_rates()
+
+# Plot rate of change (delta)
+batch.plot_deltas()
+```
+
+## Saving Results
+
+Save results to HDF5 format:
+
+```python
+batch.save_results_in_hdf5()
+```
+
+This creates `results.h5` containing:
+- `time`: Time array
+- `concentrations`: All species concentrations
+- `estimated_rates`: Calculated rates
+- `parameters`: Constants used
+
+## Advanced: Henry Equilibrium
+
+For gas-liquid partitioning:
+
+```python
+# Add both aqueous and gas species
+batch.add_species(name='CO2_aq', init_conc=0.001)
+batch.add_species(name='CO2_gas', init_conc=0)
+
+# Set Henry equilibrium (Hcc is dimensionless Henry constant)
+batch.henry_equilibrium(aq='CO2_aq', gas='CO2_gas', Hcc=0.83)
+```
+
+## Advanced: Acid-Base Equilibrium
+
+For pH-dependent systems:
+
+```python
+# Add carbonate species
+batch.add_species(name='H2CO3', init_conc=0.001)
+batch.add_species(name='HCO3', init_conc=0.01)
+batch.add_species(name='CO3', init_conc=0.0001)
+
+# Define acid with pKa values
+batch.add_acid(
+    species=['H2CO3', 'HCO3', 'CO3'],
+    pKa=[6.35, 10.33],
+    charge=0  # Charge of fully protonated form
+)
+
+# Add non-dissociating ion
+batch.add_species(name='Na', init_conc=0.1)
+batch.add_ion(name='Na', charge=1)
+
+# Create the acid-base system
+batch.create_acid_base_system()
+```
+
+After `create_acid_base_system()`, a `pH` species is automatically created and tracked.
+
+## Complete Example
+
+```python
+from porousmedialab.batch import Batch
+
+# Create batch reactor
+batch = Batch(tend=40, dt=1)
+
+# Add species
+batch.add_species(name='OM', init_conc=0.01)    # Organic matter
+batch.add_species(name='O2', init_conc=0.2)     # Oxygen
+batch.add_species(name='CO2', init_conc=0)      # Carbon dioxide
+
+# Set constants
+batch.constants['k'] = 0.5      # Degradation rate
+batch.constants['Km'] = 0.02    # Half-saturation
+
+# Define rate
+batch.rates['R_deg'] = 'k * OM * O2 / (Km + O2)'
+
+# Set time derivatives
+batch.dcdt['OM'] = '-R_deg'
+batch.dcdt['O2'] = '-R_deg'
+batch.dcdt['CO2'] = 'R_deg'
+
+# Run and plot
+batch.solve()
+batch.plot_profiles()
+```
+
+## Tips
+
+1. **Timestep selection**: Start with a larger `dt` and decrease if you see numerical instabilities
+2. **Units**: Be consistent with units throughout (e.g., mol/L for concentrations, 1/day for rates)
+3. **Debugging**: Use `verbose=True` to see simulation progress
+4. **Rate expressions**: Ensure all variables in rate expressions are defined as species or constants