Context. Magnetic fields play a pivotal role in the formation and evolution of low-mass stars, but the dynamo mechanisms generating these fields are poorly understood. Measuring cool star magnetism is a complicated task because of the complexity of cool star spectra and the subtle signatures of magnetic fields. Aims. Based on detailed spectral synthesis, we carry out quantitative measurements of the strength and complexity of surface magnetic fields in the four well-known M dwarfs GJ 388, GJ 729, GJ 285, and GJ 406 that populate the mass regime around the boundary between partially and fully convective stars. Very high-resolution (R = 100 000), high signal-to-noise (up to 400), near-infrared Stokes I spectra were obtained with CRIRES at ESO's Very Large Telescope covering regions of the FeH Wing-Ford transitions at 1 μm and Na i lines at 2.2 μm.Methods. A modified version of the Molecular Zeeman Library (MZL) was used to compute Landé g-factors for FeH lines. We determined the distribution of magnetic fields by magnetic spectral synthesis performed with the Synmast code. We tested two different magnetic geometries to probe the influence of field orientation effects. Results. Our analysis confirms that FeH lines are excellent indicators of surface magnetic fields in low-mass stars of type M, particularly in comparison to profiles of Na i lines that are heavily affected by water lines and that suffer problems with continuum normalization. The field distributions in all four stars are characterized by three distinct groups of field components, and the data are consistent neither with a smooth distribution of different field strengths nor with one average field strength covering the full star. We find evidence of a subtle difference in the field distribution of GJ 285 compared to the other three targets. GJ 285 also has the highest average field of 3.5 kG and the strongest maximum field component of 7-7.5 kG. The maximum local field strengths in our sample seem to be correlated with rotation rate. While the average field strength is saturated, the maximum local field strengths in our sample show no evidence of saturation. Conclusions. We find no difference between the field distributions of partially and fully convective stars. The one star with evidence of field distribution different from the other three is the most active star (i.e. with X-ray luminosity and mean surface magnetic field) rotating relatively fast. A possible explanation is that rotation determines the distribution of surface magnetic fields, and that local field strengths grow with rotation even in stars in which the average field is already saturated.