diff --git a/doc/source/pre/whatsnew/1.4.rst b/doc/source/pre/whatsnew/1.4.rst
index 4bb9058e9ca04c531e666e0952e6d28d4277811f..193b4056f414c8a8c3482140b1b4eda3f9f2a7c1 100644
--- a/doc/source/pre/whatsnew/1.4.rst
+++ b/doc/source/pre/whatsnew/1.4.rst
@@ -24,24 +24,24 @@ would be washed out by the presence of the peak. Now `~kwant.plotter.map`
 employs a heuristic for setting the colorscale when there are outliers,
 and will emit a warning when this is detected.
 
-System parameter names can be modified
---------------------------------------
-After the introduction of ``Builder.fill`` it has become common to construct
+System parameter substitution
+-----------------------------
+After the introduction of ``Builder.fill`` it has become possible to construct
 Kwant systems by first creating a "model" system with high symmetry and then
-filling a lower symmetry system with this model. Often, however, you want
-to use different parameter values in different parts of your system. In
+filling a lower symmetry system with this model. Often, however, one wants
+to use different parameter values in different parts of a system. In
 previous versions of Kwant this was difficult to achieve.
 
-Builders now have a method ``subs`` that makes it easy to substitute different
-names for parameters. For example if you have a Builder ``model`` that has
-a parameter ``V``, and you wish to have different values for ``V`` in your
-scattering region and leads you could do the following::
+Builders now have a method ``substitute`` that makes it easy to substitute
+different names for parameters. For example if a builder ``model``
+has a parameter ``V``, and one wishes to have different values for ``V`` in
+the scattering region and leads, one could do the following::
 
    syst = kwant.Builder()
-   syst.fill(model.subs(V='V_dot', ...))
+   syst.fill(model.substitute(V='V_dot', ...))
 
    lead = kwant.Builder()
-   lead.fill(model.subs(V='V_lead'), ...)
+   lead.fill(model.substitute(V='V_lead'), ...)
 
    syst.attach_lead(lead)
    fsyst = syst.finalized()
diff --git a/kwant/builder.py b/kwant/builder.py
index 5039f92275375589de19c6393e0cabf011c75ab8..394a701aa69efd7c9b276653ef6a132b2d484906 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1340,7 +1340,7 @@ class Builder:
         self.update(other)
         return self
 
-    def subs(self, **subs):
+    def substitute(self, **subs):
         """Return a copy of this Builder with modified parameter names.
         """
         # Get value *functions* only
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 915731e48387f016fefe34418fc9cf58416d07c7..9204ef65e93ec2212eb16f5b7b78d145e660a01a 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -1316,30 +1316,30 @@ def test_subs():
 
     syst = make_system()
     # substituting a paramter that doesn't exist produces a warning
-    warns(RuntimeWarning, syst.subs, fakeparam='yes')
+    warns(RuntimeWarning, syst.substitute, fakeparam='yes')
     # name clash in value functions
-    raises(ValueError, syst.subs, b='a')
-    raises(ValueError, syst.subs, b='c')
-    raises(ValueError, syst.subs, a='site')
-    raises(ValueError, syst.subs, c='sitea')
-    # cannot call 'subs' on systems with attached leads, because
+    raises(ValueError, syst.substitute, b='a')
+    raises(ValueError, syst.substitute, b='c')
+    raises(ValueError, syst.substitute, a='site')
+    raises(ValueError, syst.substitute, c='sitea')
+    # cannot call 'substitute' on systems with attached leads, because
     # it is not clear whether the substitutions should propagate
     # into the leads too.
     syst = make_system()
     lead = make_system(kwant.TranslationalSymmetry((-1,)), n=1)
     syst.attach_lead(lead)
-    raises(ValueError, syst.subs, a='d')
+    raises(ValueError, syst.substitute, a='d')
 
     # test basic substitutions
     syst = make_system()
     expected = hamiltonian(syst, a=1, b=2, c=3)
     # 1 level of substitutions
-    sub_syst = syst.subs(a='d', b='e')
+    sub_syst = syst.substitute(a='d', b='e')
     assert np.allclose(hamiltonian(sub_syst, d=1, e=2, c=3), expected)
     # 2 levels of substitution
-    sub_sub_syst = sub_syst.subs(d='g', c='h')
+    sub_sub_syst = sub_syst.substitute(d='g', c='h')
     assert np.allclose(hamiltonian(sub_sub_syst, g=1, e=2, h=3), expected)
     # very confusing but technically valid. 'a' does not appear in 'hopping',
     # so the signature of 'onsite' is valid.
-    sub_syst = syst.subs(a='sitea')
+    sub_syst = syst.substitute(a='sitea')
     assert np.allclose(hamiltonian(sub_syst, sitea=1, b=2, c=3), expected)