def METHOD_NAME(refinement, **model_options):
    # Domain
    lx = 10
    ly = 1
    nx = 40*refinement
    ny = 4
    mesh2d = RectangleMesh(nx, ny, lx, ly)
    depth = 40.0
    # Time interval
    dt = 0.1/refinement
    t_end = 1.0
    t_export = 0.1
    model_options['simulation_end_time'] = t_end
    # Bathymetry
    P1_2d = get_functionspace(mesh2d, 'CG', 1)
    bathy_2d = Function(P1_2d, name='Bathymetry')
    bathy_2d.assign(depth)
    # Diffusive flux BC to impose
    diff_flux = 0.2
    # Get truncated Fourier series solution
    fsol = fourier_series_solution(mesh2d, lx, diff_flux, **model_options)
    # Solver
    solver_obj = solver2d.FlowSolver2d(mesh2d, bathy_2d)
    options = solver_obj.options
    options.output_directory = 'outputs'
    options.no_exports = True
    options.timestep = dt
    options.simulation_export_time = t_export
    nu = model_options.pop('horizontal_diffusivity')
    options.add_tracer_2d('tracer_2d', 'Depth averaged tracer', 'Tracer2d',
                          diffusivity=nu)
    options.tracer_only = True
    options.horizontal_diffusivity_scale = nu
    options.horizontal_velocity_scale = Constant(0.0)
    options.fields_to_export = ['tracer_2d']
    options.update(model_options)
    options.use_limiter_for_tracers = options.tracer_element_family == 'dg'
    options.use_supg_tracer = options.tracer_element_family == 'cg'
    options.simulation_end_time = t_end - 0.5*dt
    # Boundary conditions
    solver_obj.bnd_functions['tracer_2d'] = {1: {'diff_flux': diff_flux*nu}}
    # NOTE: Zero diff_flux boundaries are enforced elsewhere by default
    # Run model
    solver_obj.assign_initial_conditions()
    solver_obj.iterate()
    sol = solver_obj.fields.tracer_2d
    if not options.no_exports:
        File('outputs/finite_element_solution.pvd').write(sol)
        File('outputs/fourier_series_solution.pvd').write(fsol)
    return errornorm(sol, fsol)