General Comments

This paper presents valuable contributions to the field of plastic modeling. Specifically, it helps address specific knowledge gaps in the vertical distribution of microplastics by focusing on biofouling, though it also includes other vertical mechanisms. This works goes even further to investigate global and seasonal variations that can impact biofilm growth vertical transport. Overall, the paper is clear and well-written, and I only recommend minor revisions before publication.

Specific Comments

The integration time step is stated to be 60 seconds around line 95. Did you examine the sensitivity of the model to this choice for time step? Is there some other justification?

There is discussion on initial density around line 100 to justify the model choice of 920 kg/m^3, but this only mentions comparison to particles with lower density (down to 30 kg/m^3). Do you know what would happen for particles closer to the density of water but still positively buoyant such as HDPE?

Around line 180 “We define the biofilm density as 1170, although the use of a denser biofilm does not change our results.” What about less dense than 1170 kg/m^3? Why would very heavy diatoms be floating at the surface? Is there a justification for that?

Figures 4 and 5: It might be good in the caption to explain more about the white bar removed section from around 400 m to 2000m. What would be going on in this region? A continuation of the behavior from 0 to 400 m?

Though Section 3.4 does a good job addressing the model assumptions, I think it would help to address a few other assumptions. This model relies on the assumption of defouling, which has only been observed in one study, and oscillations have never been experimentally observed. Additionally, this model assumes any biofilm attachment will be significantly denser than the water density. It may be worth addressing this assumption also.

Technical Corrections

Line 5: Phrase “for the physics” is awkward

Line 10: “when the processes affecting the settling velocity of the particle and the motion of the ocean are in equilibrium” it is not clear what this line means

Line 225: “our simulations only include vertical motion (advection and mixing) in order to isolate localised biological and physical effects on vertical particle displacement” this is confusing

Line 235: “This seems contrasting” is awkward, maybe “this is in contrast with the findings…”

Line 305: “hence Ggrow is the dominant down to the MLD” I think a word is missing here, dominant term?

The manuscript presents a modeling study of the vertical distribution and trajectory of microplastics at three oceanic regions characterized by contrasted biological and physical properties. The study builds on the work in Lobelle et al. (2021), including improved parameterizations for vertical mixing and biofilm loss terms. The study is a good contribution to gain insight into the mechanisms driving the vertical transport of microplastics in the ocean. The manuscript is well written, with a clear structure. Therefore, I support the publication of this work.

Some minor comments are given below:

The study is based on several model parameterizations and assumptions. The most important assumptions are discussed in Section 3.4 together with future potential developments. This section is quite pertinent, however, the model performance strongly depends on the calibration of a large number of parameters, a fact that may be further discussed or clarified :

• Parameters of mixing formulation: e.g. roughness scale, wave age, etc.
• Parameters of biofouling formulations: e.g. collision rate, growth rate, etc
• Would the combined use of a different plastic density and biofilm density change the results?
• How is the calibration/validation of NEMO-MEDUSA? Is it sensible to different parameterizations as well?

Kooi’s model seems to depend on water physical parameters such as temperature, salinity, viscosity. Did this work consider the seasonal and spatial variability of these parameters or are they just affecting the biofouling parameters in NEMO-MEDUSA? If so, how do they affect the spatial differences of the vertical distribution of particles between the three regions?

Unless I am mistaken, the discussion of results from Appendix D is quite light. Given that the oscillatory behavior of microplastics has not been observed yet, it may be interesting to elaborate on this scenario that may also represent a "realistic" situation.

Line 25: A new study published in Science (Weiss et al., 2021) reformulates the calculation of plastic fluxes and shows that there may not be this huge amount of “missing plastic”.

Weiss, L., Ludwig, W., Heussner, S., Canals, M., Ghiglione, J.F., Estournel, C., Constant, M., Kerhervé, P. (2021). The missing ocean plastic sink: Gone with the rivers, Science, 373 (6550), 107-111

This is a very interesting paper building on previous work by some of the authors, notably, Kooi et al. (2017) and Lobelle et al. (2021). The authors investigate the dynamics of submerged plastic particles in response to vertical advection, diffusion and biofouling using for that a well-established Lagrangean model of marine plastic transport. Main results of this study are that plastic particles still undergo oscillations, similar to those discussed by Kooi et al. (2017) and Lobelle et al. (2021), when the effects of vertical advection and diffusion are taken into account. Subsurface maximum concentrations of plastics can be created in low turbulence areas when downward settling is approximately balance by ambient upwelling. Strong mixing can draw plastic particles down to hundreds of metres below the surface, for example in the Southern Ocean.

The paper is well written, well referenced for the most part, and the arguments put forward to explain the model results are generally rigorous and persuasive. I have nevertheless quite a few comments and remarks to make, and I have some objections regarding aspects of the methodology and interpretation of the results.  This makes me recommend major revisions. In the following, my comments are organised by line number.

Ll. 41-45. Vertical mixing within the mixed layer is driven both by wind stirring and buoyancy losses. While the Kulkuka et al. (2012) parameterisation accounts only for wind mixing, the mixing that takes place in the NEMO-MEDUSA and that accounts in part, for the vertical distribution of phytoplankton and primary production the authors use in their simulations, takes into account both wind stirring and instabilities driven by buoyancy losses. Similarly, vertical mixing of plastic particles in the mixed layer is caused by both turbulent mechanisms, not just wind-induced mixing. In a similarly vein, mixing in the ocean interior stems from many causes, one of them being the dissipation of tidal energy on the shelves (this mixing will often overlap with mixing originating from surface buoyancy and momentum fluxes) and in the deep ocean over rough topography. But dissipation of tidal energy is not the sole source of turbulent kinetic energy in the deep ocean and, away from boundaries and rough bathymetry, not even the most important. So, it is not entirely appropriate to refer to ocean interior mixing as tidal mixing. Finally, Mountford and Morales Maqueda (2019, 2020) have also demonstrated the importance of diapycnal mixing through the water column in controlling the penetration of buoyant microplastics in the water column. Please acknowledge.

Ll. 52-54. That these oscillations may exist is entirely plausible but, to my knowledge, they have never been identified in the field. One of the beauties of modelling is that it allows us to explore and predict phenomena that have not yet been observed. Given the growing amount of literature on the subject of these plastic oscillations, it is becoming urgent to come up with ideas on how to gather observational evidence of the process to ascertain whether it is actually taking place in the real world.

Ll 95-96. The validity of Stokes’ law does depend on the radius of the particles but on the Reynolds number. Please, correct and provide the correct reason, which much be that buoyant particles with a larger radius will shift the motion regime to one of high Reynolds numbers for with Stokes’ law needs to be corrected.

Ll 100-105. Are you sure that the use of particles as light as  30 kg m^-3 does not push the system to terminal Reynolds number much smaller than 1, outside the permissible Stokes regime?

Ll 115-117. One-dimensionality is a dangerous assumption here. Water is nearly incompressible, and incompressible fluid motion is never one-dimensional if the divergence of the flow in any one direction is non-zero. This must be, in general, the case in your simulations. For example, an annual mean upward vertical velocity of 10^-9 m s^-1 (Fig. A1) at the base of an equatorial column of water extending over a 1^oÎ1^o area and 1000 m deep would cause a horizontal divergent current of about 1 cm s^-1, which means that any particle of plastic originally within the 1^oÎ1^o would have left the area in about one hundred days.

Ll 118-120. This is a rather bizarre statement. It is not clear to what is meant by “very low-resolution of turbulence”. In ocean models, the resolution of diapycnal turbulence is as high or as low as that of any other vertical or quasi-vertical process. If you discard the turbulence calculated by the NEMO-MEDUSA model because they are allegedly of too low resolution, it is not clear to me why would you think acceptable to use the vertical velocities produced by the same model and which are computed on the same low resolution, especially since the dynamics of turbulence and vertical advection are intimately related. The formulation of turbulent vertical diffusivities in NEMO-MEDUSA is based in the TKE scheme of Gaspar et al. (1990) which is one of the best physically founded and conceptually appealing formulation of these processes for ocean modelling. Like the KPP formulation, the TKE scheme is capable of representing non-uniform turbulent profiles in the simulated mixed layer as long as the said mixed layer is resolved well enough by the model. This is not to denigrate the author’s choice of the KPP parameterisation to calculate diffusivities in their model. I presume this is because they do not have access to the diffusivities produced by the NEMO-MEDUSA model. This is fine as a justification for your approach. But then, it suffices to state so. There is no need to advance questionable arguments about the resolution of turbulence in NEMO-MEDUSA.

L 125 u_star is called the friction, not frictional, velocity (see also L 130). It does not correspond to any velocity of the seawater at the surface or anywhere else. It is simply a way of expressing a stress, in this case the surface wind stress, in terms of a velocity.

L 127. Ignoring the impact of Langmuir enhanced mixing in the Southern Ocean might be problematic but perhaps acceptable here for simplicity.

Ll 133-135. This statement is also rather obscure. The non-locality of the KPP scheme has nothing to do with its (inexistent) dependence on the wind stress at different locations. It refers to the fact that, in cases of unstable forcing, turbulent fluxes are decomposed into a local term that is proportional to the vertical gradient of the property being fluxed and a non-local term that conveys information of the surface property flux down through the boundary layer. There is, however, nothing wrong with neglecting this nonlocal term, again, for simplicity.

L 143. It is Delta t, not delta t

Ll 146-148. Ocean interior mixing is driven by tidal processes as well as by processes other than tides. Instead of an extraneous set of diffusivities, it might have been desirable to use the diffusivity output from NEMO-MEDUSA, which includes all relevant sources of turbulent forcing and, very importantly, is physically compatible with the vertical velocities you are using. I appreciate this might not be possible, though. In any case, could the authors please show the vertical diffusivity profiles (in annual mean, say) they are using in the simulations for the three sites?

Ll 164-166 Please expand. In which way is the assumption of diatom-dominated biofilms better than one in which both diatoms and non-diatoms are combined?

Ll 236-238. As indicated above, other authors have, however, demonstrated the importance of advection and diffusion processes, both diapycnal (~vertical) and isopycnal (~horizontal), for the distribution of microplastics in the ocean.

L 265. I find this phrase awkward, could you please rephrase?

L 288. I am not sure this is the most accurate way of describing what is happening. The particles in the mixed layer are being strongly mixed (passively), rather than moving “passively with the flow”.

Figure 4 caption. What, pray tell, is Wmixing? Is this given by (3). If so, use a consistent notation.

1. 312 A potential reason? Please, look at the data of primary production you are using and confirm or disprove your hypothesis. There is no reason here to advance a guess when the actual explanation is at hand.

L 325 Equation (5) does not readily show that larger particles have larger settling velocities than smaller ones, although your statement is correct as long as the densities of the particles are the same.

L 343 heterotrophs.

1. 409 Lresp is repeated twice.