Part 2: Solving a non-linear initial-value problem with the implicit Euler scheme

A municipality plans to build a large battery energy storage facility and is looking for a suitable location. The potential location is close (1 km) to a well that is used to abstract drinking water from an aquifer with a constant pumping rate \(Q\). In the event of a fire or accident in the battery energy storage facility, contaminants could infiltrate into the subsurface. Therefore, the municipality wants to assess the risk of contamination in case a leakage occurs. One requirement is that the facility cannot be built in the groundwater protection zone, which is delineated by the area in which a water parcel takes less than 25 years to reach the well.

Assuming radial flow towards the well in a uniform aquifer, the distance \(r\) of a groundwater particle from the well is governed by

\[\frac{dr}{dt} = - \frac{Q}{2\pi r n b} = -\alpha\frac{1}{r}\,,\]

where \(n\) is porosity, \(b\) is the aquifer thickness, and \(\alpha = \frac{Q}{2\pi nb}\).

To find out how long it takes for a particle starting at a distance \(r_0\) to reach the well, you need to solve the ODE until the distance \(r\) equals the well radius, \(r_{\mathrm{well}}\).

import numpy as np
import matplotlib.pyplot as plt

Q = 3000 # m³/d
b = 20 # m
n = 0.3

alpha = Q / (2 * np.pi * b * n)


r0 = 1000 # m
r_well = 0.2 # m

Task 2.1:

Formulate the implicit Euler scheme for the ODE.

Solution 2.1:

\[r_{i+1} = r_{i} + \Delta t (-\alpha) \frac{1}{r_{i+1}}\]

For general nonlinear ODEs, the implicit Euler update cannot usually be solved directly. Instead, we use an iterative scheme – Newton-Raphson – to solve for \(r_{i+1}\).¹

In the Newton-Raphson scheme, the expression for the new value, \(r_{i+1}^{j+1}\), is given by:

\[r_{i+1}^{j+1} = r_{i+1}^{j} - \frac{q(r_{i+1}^{j})}{q^\prime(r_{i+1}^{j})}\]

Here \(q\) is obtained by rearranging the implicit Euler expression so that one side of the equation is zero.

¹ Note that the equation in this specific example is quadratic (after multiplying both sides with \(r_{i+1}\)), and thus, could be solved directly, even though it is non-linear. For the sake of practicing, we ask you to use the Newton-Raphson method nevertheless.

Task 2.2:

Formulate an expression for \(q\) by rearranging the equation for the implicit Euler scheme from task 2.1.

Solution 2.2

\[q(r_{i+1}) = r_{i+1} - r_{i} - \Delta t (-\alpha) \frac{1}{r_{i+1}}\]

For the Newton Raphson scheme we also need its derivative \(q'\):

$$q’(r_{i+i}) = \frac{\mathrm{d}q}{\mathrm{d}r_{i+1}} = 1 - \Delta t \alpha\frac{1}{r_{i+1}^2}

Task 2.3:

Complete the code cell below to implement an implicit Euler scheme with Newton Raphson iterations for the particle-tracking ODE.

Note: the version of the notebook distributed in class contained some errors in the code template that were fixed now.

• r_next was not updated in the inner loop

• The convergence test checked if r_update did not change anymore instead of checking that q becomes small enough.

# Solve the IVP with an implicit Euler scheme, using the Newton-Raphson method

# Initialize variables
r_next = r0
r_previous = r0
t_previous = 0
dt = 50 # days

times = [t_previous]
r = [r0]

# Set parameters for the Newton-Raphson method and time iteration
tolerance = 1e-6
max_times = 1e4
max_iter = 50

# Loop over time
for iteration in range(int(max_times)):
    t_next = t_previous + dt
    times.append(t_next)

    # Compute q and q_prime once for the current guess r_next
    q = r_next - r_previous - dt * (-alpha / r_next)
    q_prime = 1 - dt * alpha / r_next ** 2
    iter_count = 0

    # Loop for Newton-Raphson iterations
    while iter_count < max_iter:
        # Update the guess using Newton-Raphson formula
        r_update = r_next - q / q_prime

        # Compute residual and derivative at the updated guess
        q_update = r_update - r_previous - dt * (-alpha / r_update)
        q_prime_update = 1 - dt * alpha / r_update ** 2

        # Check for convergence on the residual (root)
        if abs(q_update) < tolerance:
            r_next = r_update
            break

        # Promote the updated values for the next iteration (avoid recomputing q)
        r_next = r_update
        q = q_update
        q_prime = q_prime_update
        iter_count += 1


    if iter_count >= max_iter:
        print(f"Warning: Newton did not converge at t={t_next} (iter={iter_count})")
    
    # Update the values for the next time step using the converged r_next
    r_previous = r_next
    r.append(r_next)
    t_previous = t_next

    # Stop the simulation when the particle reaches the well (r <= 0.2 m)
    if r_next <= r_well:
        break

Warning: Newton did not converge at t=6200 (iter=50)

Solution 2.3

Note: The implementation above is intentionally kept simple to illustrate the core algorithm. In practical applications, additional safeguards and refinements are often required to improve robustness and efficiency:

• The derivative \(q^\prime\) may become zero when \(r_{i+1}^j = \sqrt{\Delta t \alpha}\). In this case, the Newton update is undefined, causing the iteration to fail. A common safeguard is to reduce the time step when \(|q'|\) becomes too small, or to switch temporarily to a more robust root-finding strategy.

• As the particle approaches the well, the velocity increases because of the \(\frac{1}{r}\)​ term. With a fixed time step, a large \(\Delta t\) may cause the simulation to overshoot the point at which the particle reaches the well, while a very small \(\Delta t\) can make the computation unnecessarily expensive when the particle is still far away. Adaptive time stepping, in which the time step is adjusted based on the behaviour of the solution, can improve both accuracy and efficiency.

• The Newton iteration uses only an absolute convergence tolerance. In practice, it is often preferable to combine absolute and relative tolerances. A relative tolerance accounts for the scale of the solution, preventing convergence criteria from becoming either too strict for large values or too loose for small values.

Task 2.4:

Run the code cell below to plot the numerical solution against the analytical solution

years = np.array(times) / 365.25
plt.plot(years, r, label="implicit Euler scheme")
plt.xlabel("time [years]")
plt.ylabel("distance from the well [m]")
r_analytical = np.sqrt(-2 * alpha * np.array(times) + r0**2)
plt.plot(years, r_analytical, label="analytical solution")
plt.legend()
plt.ylim(bottom=0, top =r0)

(0.0, 1000.0)

Task 2.5:

Run the code cell below to compare the time it takes until the particle reaches the well between the numerical and analytical solution.

t_end_analytical = (r0**2 - r_well**2) / 2 / alpha
t_end_numerical = times[-1]

print("Time needed to reach the well (r=r_well):\n")
print(f"Analytical solution:\t{t_end_analytical / 365.25:.2f} years")
print(f"Numerical solution:\t{t_end_numerical / 365.25:.2f} years")

Time needed to reach the well (r=r_well):

Analytical solution:	17.20 years
Numerical solution:	16.97 years

Task 2.6:

Based on this estimate, would you recommend approving the proposed location? Explain your reasoning.

Solution 2.6

Since the travel time to the well is about 17 years, the battery energy storage facility falls into the protection zone, which is based on the 25-year contour line.

Based on this first estimate, the proposed location cannot be approved. Even though the model is strongly simplified, it is likely not worth setting up a more complex model the model that accounts for, e.g., dispersion, sorption, degradation, regional groundwater flow and aquifer heterogeneity.

By Ronald Brinkgreve, Anna Störiko, Delft University of Technology. CC BY 4.0, more info on the Credits page of Workbook.